diff --git a/src/meshmodifiers/InterfaceFromSideset.C b/src/meshmodifiers/InterfaceFromSideset.C index c2d4217c..52c4bd26 100644 --- a/src/meshmodifiers/InterfaceFromSideset.C +++ b/src/meshmodifiers/InterfaceFromSideset.C @@ -69,9 +69,12 @@ InterfaceFromSideset::modify() _mesh_ptr->buildBndElemList(); - // Identify node IDs on given boundaries + // Identify node IDs on given boundaries and all sidesets std::set boundary_node_ids; + std::set allsidesets_node_ids; + std::map> sidesets_node_ids; for (auto it = _mesh_ptr->bndElemsBegin(); it != _mesh_ptr->bndElemsEnd(); ++it) + { if (boundary_ids.count((*it)->_bnd_id) > 0) { Elem * elem = (*it)->_elem; @@ -80,6 +83,22 @@ InterfaceFromSideset::modify() for (unsigned int n = 0; n < nodes_on_side.size(); ++n) boundary_node_ids.insert(elem->node_id(nodes_on_side[n])); } + for (auto & sideset_name : sideset_names) + { + auto sidedset_id = _mesh_ptr->getBoundaryID(sideset_name); + if ((*it)->_bnd_id == sidedset_id) + { + Elem * elem = (*it)->_elem; + auto s = (*it)->_side; + std::vector nodes_on_side = elem->nodes_on_side(s); + for (unsigned int n = 0; n < nodes_on_side.size(); ++n) + { + sidesets_node_ids[sideset_name].insert(elem->node_id(nodes_on_side[n])); + allsidesets_node_ids.insert(elem->node_id(nodes_on_side[n])); + } + } + } + } // identify nodes "not at the edge" for each sideset std::map> nb_neighbors_on_sideset; @@ -118,17 +137,12 @@ InterfaceFromSideset::modify() std::set treated_node_ids; _new_boundary_sides_map.clear(); - // identify all nodes on this sideset - std::set sideset_node_ids; - for (auto it = _mesh_ptr->bndElemsBegin(); it != _mesh_ptr->bndElemsEnd(); ++it) - if ((*it)->_bnd_id == sidedset_id) - { - Elem * elem = (*it)->_elem; - auto s = (*it)->_side; - std::vector nodes_on_side = elem->nodes_on_side(s); - for (unsigned int n = 0; n < nodes_on_side.size(); ++n) - sideset_node_ids.insert(elem->node_id(nodes_on_side[n])); - } + // get a set of nodes on all sidesets but this sideset + std::set othersidesets_node_ids; + for (auto & sideset_name2 : sideset_names) + if (sideset_name2 != sideset_name) + othersidesets_node_ids.insert(sidesets_node_ids[sideset_name2].begin(), + sidesets_node_ids[sideset_name2].end()); // loop on all elements involved in that sideset for (auto it = _mesh_ptr->bndElemsBegin(); it != _mesh_ptr->bndElemsEnd(); ++it) @@ -147,16 +161,19 @@ InterfaceFromSideset::modify() const bool is_node_on_boundary = boundary_node_ids.find(current_node_id) != boundary_node_ids.end(); const bool has_two_neighbours = nb_neighbors_on_sideset[sideset_name][current_node_id]>2; - //printf(" Node %d has %d neighbors on sideset\n", current_node_id, nb_neighbors_on_sideset[sideset_name][current_node_id]); + const bool is_node_on_other_sidesets = + othersidesets_node_ids.find(current_node_id) != othersidesets_node_ids.end(); bool do_duplicate = false; if (is_node_on_boundary) do_duplicate = true; else if (has_two_neighbours) do_duplicate = true; + else if (is_node_on_other_sidesets) + do_duplicate = true; if (do_duplicate) { // node to be duplicated - printf(" Node %d needs to be duplicated\n", current_node_id); + printf(" Node %d needs to be duplicated", current_node_id); const Node * current_node = mesh.node_ptr(current_node_id); std::set elem_ids_on_that_side; std::set elem_ids_on_other_side; @@ -166,10 +183,12 @@ InterfaceFromSideset::modify() new_node = Node::build(*current_node, mesh.n_nodes()).release(); new_node->processor_id() = current_node->processor_id(); mesh.add_node(new_node); + printf(" -> node %d\n", new_node->id()); if (is_node_on_boundary) boundary_node_ids.insert(new_node->id()); treated_node_ids.insert(new_node->id()); // to avoid looping on new node - sideset_node_ids.insert(new_node->id()); + sidesets_node_ids[sideset_name].insert(new_node->id()); + allsidesets_node_ids.insert(new_node->id()); std::map>::const_iterator node_to_elem_pair = node_to_elem_map.find(current_node_id); if (node_to_elem_pair==node_to_elem_map.end()) @@ -177,8 +196,13 @@ InterfaceFromSideset::modify() new_node_to_elem_map[new_node->id()] = node_to_elem_pair->second; for (auto & sideset_name2 : sideset_names) if (sideset_name2 != sideset_name) + { nb_neighbors_on_sideset[sideset_name2][new_node->id()] = nb_neighbors_on_sideset[sideset_name2][current_node_id]; + if (sidesets_node_ids[sideset_name2].find(current_node_id) + != sidesets_node_ids[sideset_name2].end()) + sidesets_node_ids[sideset_name2].insert(new_node->id()); + } // Add sideset/boundary info to the new node std::vector node_boundary_ids = @@ -198,7 +222,7 @@ InterfaceFromSideset::modify() // Find which side is it from our segment bool is_on_that_side = isElementOnThatSideOfSegment( mesh, current_elem, current_node_id, other_node_id, - sideset_node_ids, normal_vec); + sidesets_node_ids[sideset_name], normal_vec); if (is_on_that_side) { elem_ids_on_that_side.insert(elem_id); @@ -269,12 +293,6 @@ InterfaceFromSideset::isElementOnThatSideOfSegment( } } -// printf("isElementOnThatSideOfSegment on element %d ?\n",elem->id()); -// if (is_node1_smaller_than_node2) -// printf(" Ordered segment: nodes %d - %d\n", node_id1,node_id2); -// else -// printf(" Ordered segment: nodes %d - %d\n", node_id2,node_id1); - bool result = true; // element is on "that" (one) side of segment for (MooseIndex(elem->n_vertices()) i_v = 0; i_v < elem->n_vertices(); ++i_v) { @@ -290,8 +308,7 @@ InterfaceFromSideset::isElementOnThatSideOfSegment( else is_node_on_that_side = isNodeOnThatSideOfSegment( test_node, node2, node1, normal_vec); - //printf(" isElement %d OnThatSideOfSegment? %d\n", elem->id(), is_node_on_that_side); - return is_node_on_that_side; // one node on one side is enough to decide + return is_node_on_that_side; // one node on one side is enough to decide } } mooseError("No node evaluated in isElementOnThatSideOfSegment"); @@ -323,7 +340,6 @@ InterfaceFromSideset::isNodeOnThatSideOfSegment( bool result = false; if (test > 0) result = true; - //printf(" isNode %d OnThatSideOfSegment %d-%d? %d\n",nodetest.id(),node1.id(),node2.id(), result); return result; } @@ -343,7 +359,6 @@ InterfaceFromSideset::getMeshNormalVector(const MeshBase & mesh) * ((nodeC)((i+2)%LIBMESH_DIM) - (nodeA)((i+2)%LIBMESH_DIM)) - ((nodeB)((i+2)%LIBMESH_DIM) - (nodeA)((i+2)%LIBMESH_DIM)) * ((nodeC)((i+1)%LIBMESH_DIM) - (nodeA)((i+1)%LIBMESH_DIM)); - //printf(" Found normal vector %f, %f, %f\n",normal_vec[0],normal_vec[1],normal_vec[2]); return normal_vec; } diff --git a/tests/meshmodifiers/fractureX.geo b/tests/meshmodifiers/fractureX.geo index fcfa842c..56d13d80 100644 --- a/tests/meshmodifiers/fractureX.geo +++ b/tests/meshmodifiers/fractureX.geo @@ -49,7 +49,7 @@ Physical Surface("blockD") = {4}; Physical Curve("ss2") = {5,6}; Physical Curve("ss3") = {7,8}; Physical Curve("bottom") = {2,3,4}; -Physical Curve("top") = {10,11}; +Physical Curve("top") = {10,11,12}; Physical Curve("right") = {9}; Physical Curve("left") = {1}; diff --git a/tests/meshmodifiers/fractureY.geo b/tests/meshmodifiers/fractureY.geo new file mode 100644 index 00000000..5f49e2c6 --- /dev/null +++ b/tests/meshmodifiers/fractureY.geo @@ -0,0 +1,52 @@ +// Generate simple X across 2D box + +lc = 0.3; +lc2 = 0.3; + +Point(1) = {0, 0, 0, lc}; +Point(2) = {1, 0, 0, lc}; +Point(3) = {1, 1, 0, lc}; +Point(4) = {0, 1, 0, lc}; + +Point(5) = {0.5, 0.5, 0, lc2}; + +Point(6) = {0.3, 0, 0, lc2}; +Point(7) = {0.7, 1, 0, lc2}; +Point(9) = {0.25, 1, 0, lc2}; + +Line(1) = {4, 1}; +Line(2) = {1, 6}; +Line(3) = {6, 2}; +Line(9) = {2, 3}; +Line(10) = {3, 7}; +Line(11) = {7, 9}; +Line(12) = {9, 4}; + +Line(5) = {6, 5}; +Line(6) = {5, 7}; +Line(7) = {9, 5}; + + +Curve Loop(1) = {1, 2, 5, -7, 12}; +Plane Surface(1) = {1}; +Curve Loop(2) = {3, 9, 10, -6, -5}; +Plane Surface(2) = {2}; +Curve Loop(3) = {6, 11, 7}; +Plane Surface(3) = {3}; + + +Physical Surface("blockA") = {1}; +Physical Surface("blockB") = {2}; +Physical Surface("blockC") = {3}; + + +Physical Curve("ss2") = {5,6}; +Physical Curve("ss3") = {7}; +Physical Curve("bottom") = {2,3}; +Physical Curve("top") = {10,11,12}; +Physical Curve("right") = {9}; +Physical Curve("left") = {1}; + + + + diff --git a/tests/meshmodifiers/fractureY.msh b/tests/meshmodifiers/fractureY.msh new file mode 100644 index 00000000..aa24adb0 --- /dev/null +++ b/tests/meshmodifiers/fractureY.msh @@ -0,0 +1,115 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$PhysicalNames +9 +1 4 "ss2" +1 5 "ss3" +1 6 "bottom" +1 7 "top" +1 8 "right" +1 9 "left" +2 1 "blockA" +2 2 "blockB" +2 3 "blockC" +$EndPhysicalNames +$Nodes +30 +1 0 0 0 +2 1 0 0 +3 1 1 0 +4 0 1 0 +5 0.5 0.5 0 +6 0.3 0 0 +7 0.7 1 0 +8 0.25 1 0 +9 0 0.7500000000003477 0 +10 0 0.5000000000020616 0 +11 0 0.2500000000010419 0 +12 0.5333333333326683 0 0 +13 0.7666666666659077 0 0 +14 0.3999999999997372 0.2499999999993431 0 +15 0.5999999999999758 0.7499999999999395 0 +16 0.375 0.75 0 +17 1 0.2499999999994109 0 +18 1 0.4999999999986918 0 +19 1 0.7499999999993401 0 +20 0.47500000000064 1 0 +21 0.2499999999998144 0.4549999999996932 0 +22 0.1999999999995926 0.184999999999782 0 +23 0.1573932185564783 0.6586067814439295 0 +24 0.1564786437112957 0.8317213562888555 0 +25 0.6678976706301279 0.3260600841690262 0 +26 0.7765233863124585 0.5737228691665684 0 +27 0.8078085088074329 0.8148696805270933 0 +28 0.8253883951420575 0.1796487987997619 0 +29 0.8641727018358398 0.3726276422386496 0 +30 0.6386572131540997 0.1511417765936262 0 +$EndNodes +$Elements +64 +1 1 2 9 1 4 9 +2 1 2 9 1 9 10 +3 1 2 9 1 10 11 +4 1 2 9 1 11 1 +5 1 2 6 2 1 6 +6 1 2 6 3 6 12 +7 1 2 6 3 12 13 +8 1 2 6 3 13 2 +9 1 2 4 5 6 14 +10 1 2 4 5 14 5 +11 1 2 4 6 5 15 +12 1 2 4 6 15 7 +13 1 2 5 7 8 16 +14 1 2 5 7 16 5 +15 1 2 8 9 2 17 +16 1 2 8 9 17 18 +17 1 2 8 9 18 19 +18 1 2 8 9 19 3 +19 1 2 7 10 3 7 +20 1 2 7 11 7 20 +21 1 2 7 11 20 8 +22 1 2 7 12 8 4 +23 2 2 1 1 10 11 21 +24 2 2 1 1 21 11 22 +25 2 2 1 1 5 16 21 +26 2 2 1 1 21 16 23 +27 2 2 1 1 1 6 22 +28 2 2 1 1 11 1 22 +29 2 2 1 1 14 5 21 +30 2 2 1 1 14 21 22 +31 2 2 1 1 16 8 24 +32 2 2 1 1 9 10 23 +33 2 2 1 1 10 21 23 +34 2 2 1 1 4 9 24 +35 2 2 1 1 8 4 24 +36 2 2 1 1 6 14 22 +37 2 2 1 1 23 16 24 +38 2 2 1 1 9 23 24 +39 2 2 2 2 18 19 26 +40 2 2 2 2 26 19 27 +41 2 2 2 2 5 14 25 +42 2 2 2 2 6 12 14 +43 2 2 2 2 25 14 30 +44 2 2 2 2 5 25 26 +45 2 2 2 2 3 7 27 +46 2 2 2 2 14 12 30 +47 2 2 2 2 15 5 26 +48 2 2 2 2 19 3 27 +49 2 2 2 2 15 26 27 +50 2 2 2 2 2 17 28 +51 2 2 2 2 13 2 28 +52 2 2 2 2 7 15 27 +53 2 2 2 2 18 26 29 +54 2 2 2 2 26 25 29 +55 2 2 2 2 17 18 29 +56 2 2 2 2 12 13 30 +57 2 2 2 2 28 17 29 +58 2 2 2 2 25 28 29 +59 2 2 2 2 13 28 30 +60 2 2 2 2 28 25 30 +61 2 2 3 3 15 7 20 +62 2 2 3 3 16 15 20 +63 2 2 3 3 16 20 8 +64 2 2 3 3 5 15 16 +$EndElements diff --git a/tests/meshmodifiers/gold/testInterfaceFromSidesetDiffusion4.e b/tests/meshmodifiers/gold/testInterfaceFromSidesetDiffusion4.e new file mode 100644 index 00000000..1629596a Binary files /dev/null and b/tests/meshmodifiers/gold/testInterfaceFromSidesetDiffusion4.e differ diff --git a/tests/meshmodifiers/gold/testInterfaceFromSidesetDiffusion5.e b/tests/meshmodifiers/gold/testInterfaceFromSidesetDiffusion5.e new file mode 100644 index 00000000..2fa3007b Binary files /dev/null and b/tests/meshmodifiers/gold/testInterfaceFromSidesetDiffusion5.e differ diff --git a/tests/meshmodifiers/testInterfaceFromSidesetDiffusion4.i b/tests/meshmodifiers/testInterfaceFromSidesetDiffusion4.i new file mode 100644 index 00000000..dbbc20a4 --- /dev/null +++ b/tests/meshmodifiers/testInterfaceFromSidesetDiffusion4.i @@ -0,0 +1,50 @@ +[Mesh] + type = FileMesh + file = fractureY.msh +[] + +[MeshModifiers] + [break] + type = InterfaceFromSideset + boundaries = 'top bottom left right' + sidesets = 'ss2 ss3' + split_interface = false + [] +[] + +[Variables] + [dymmy_var] + [] +[] + +[Kernels] + [diff] + type = Diffusion + variable = dymmy_var + [] +[] + +[BCs] + [left] + type = DirichletBC + variable = dymmy_var + boundary = 'left' + value = 1 + [] + [right] + type = DirichletBC + variable = dymmy_var + boundary = 'right' + value = 0 + [] +[] + +[Executioner] + type = Steady + solve_type = NEWTON +[] + +[Outputs] + file_base = testInterfaceFromSidesetDiffusion4 + exodus = true +[] diff --git a/tests/meshmodifiers/testInterfaceFromSidesetDiffusion5.i b/tests/meshmodifiers/testInterfaceFromSidesetDiffusion5.i new file mode 100644 index 00000000..6b493467 --- /dev/null +++ b/tests/meshmodifiers/testInterfaceFromSidesetDiffusion5.i @@ -0,0 +1,50 @@ +[Mesh] + type = FileMesh + file = fractureY.msh +[] + +[MeshModifiers] + [break] + type = InterfaceFromSideset + boundaries = 'top bottom left right' + sidesets = 'ss2 ss3' + split_interface = false + [] +[] + +[Variables] + [dymmy_var] + [] +[] + +[Kernels] + [diff] + type = Diffusion + variable = dymmy_var + [] +[] + +[BCs] + [left] + type = DirichletBC + variable = dymmy_var + boundary = 'left' + value = 0 + [] + [right] + type = DirichletBC + variable = dymmy_var + boundary = 'right' + value = 1 + [] +[] + +[Executioner] + type = Steady + solve_type = NEWTON +[] + +[Outputs] + file_base = testInterfaceFromSidesetDiffusion5 + exodus = true +[] diff --git a/tests/meshmodifiers/tests b/tests/meshmodifiers/tests index 616d9d4a..c2c6d407 100644 --- a/tests/meshmodifiers/tests +++ b/tests/meshmodifiers/tests @@ -24,4 +24,16 @@ input = 'testInterfaceFromSidesetDiffusion3.i' exodiff = 'testInterfaceFromSidesetDiffusion3.e' [../] + [./test_interfaceFromSideset_diff4] + # Testing diffusion on Y reaching sides from left to right + type = 'Exodiff' + input = 'testInterfaceFromSidesetDiffusion4.i' + exodiff = 'testInterfaceFromSidesetDiffusion4.e' + [../] + [./test_interfaceFromSideset_diff5] + # Testing diffusion on Y reaching sides from right to left + type = 'Exodiff' + input = 'testInterfaceFromSidesetDiffusion5.i' + exodiff = 'testInterfaceFromSidesetDiffusion5.e' + [../] []