diff --git a/include/meshmodifiers/InterfaceFromSideset.h b/include/meshmodifiers/InterfaceFromSideset.h index 5f7ecb88..5b53daa6 100644 --- a/include/meshmodifiers/InterfaceFromSideset.h +++ b/include/meshmodifiers/InterfaceFromSideset.h @@ -34,7 +34,6 @@ class InterfaceFromSideset : public BreakMeshByBlockBase bool isElementOnThatSideOfSegment(const MeshBase &, const Elem *, const dof_id_type, const dof_id_type, const std::set, const std::vector); std::vector getMeshNormalVector(const MeshBase &); - std::set> _neighboring_block_list; std::map, std::set>> _new_boundary_sides_map; diff --git a/meshes/fractureX.geo b/meshes/fractureX.geo new file mode 100644 index 00000000..fcfa842c --- /dev/null +++ b/meshes/fractureX.geo @@ -0,0 +1,58 @@ +// 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(8) = {0.75, 0, 0, lc2}; +Point(9) = {0.25, 1, 0, lc2}; + +Line(1) = {4, 1}; +Line(2) = {1, 6}; +Line(3) = {6, 8}; +Line(4) = {8, 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}; +Line(8) = {5, 8}; + + +Curve Loop(1) = {1, 2, 5, -7, 12}; +Plane Surface(1) = {1}; +Curve Loop(2) = {3, -8, -5}; +Plane Surface(2) = {2}; +Curve Loop(3) = {4, 9, 10, -6, 8}; +Plane Surface(3) = {3}; +Curve Loop(4) = {6, 11, 7}; +Plane Surface(4) = {4}; + + +Physical Surface("blockA") = {1}; +Physical Surface("blockB") = {2}; +Physical Surface("blockC") = {3}; +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("right") = {9}; +Physical Curve("left") = {1}; + + + + diff --git a/meshes/fractureXincluded.geo b/meshes/fractureXincluded.geo new file mode 100644 index 00000000..e0c17f38 --- /dev/null +++ b/meshes/fractureXincluded.geo @@ -0,0 +1,43 @@ +// Generate truncated X pattern across 2D box +lc = 0.3; +lc2 = 0.2; + +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.2, 0, lc2}; +Point(7) = {0.7, 0.8, 0, lc2}; +Point(8) = {0.7, 0.1, 0, lc2}; +Point(9) = {0.3, 0.9, 0, lc2}; + +Line(1) = {4, 1}; +Line(2) = {1, 2}; +Line(9) = {2, 3}; +Line(10) = {3, 4}; + +Line(5) = {6, 5}; +Line(6) = {5, 7}; +Line(7) = {9, 5}; +Line(8) = {5, 8}; + + +Curve Loop(1) = {1, 2, 9, 10}; +Plane Surface(1) = {1}; +SetFactory("OpenCASCADE"); +Curve{5, 6, 7, 8} In Surface{1}; +Physical Surface("block") = {1}; + +Physical Curve("ss2") = {5,6}; +Physical Curve("ss3") = {7,8}; +Physical Curve("bottom") = {2}; +Physical Curve("top") = {10}; +Physical Curve("right") = {9}; +Physical Curve("left") = {1}; + + + + diff --git a/meshes/fractureXtruncated.geo b/meshes/fractureXtruncated.geo new file mode 100644 index 00000000..f67fbbc0 --- /dev/null +++ b/meshes/fractureXtruncated.geo @@ -0,0 +1,46 @@ +// Generate truncated X pattern 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.6, 0.75, 0, lc2}; +Point(8) = {0.75, 0, 0, lc2}; +Point(9) = {0.25, 1, 0, lc2}; + +Line(1) = {4, 1}; +Line(2) = {1, 6}; +Line(3) = {6, 8}; +Line(4) = {8, 2}; +Line(9) = {2, 3}; +Line(10) = {3, 9}; +Line(11) = {9, 4}; + +Line(5) = {6, 5}; +Line(6) = {5, 7}; +Line(7) = {9, 5}; +Line(8) = {5, 8}; + + +Curve Loop(1) = {1, 2, 3, 4, 9, 10, 11}; +Plane Surface(1) = {1}; +SetFactory("OpenCASCADE"); +Curve{5, 6, 7, 8} In Surface{1}; +Physical Surface("block") = {1}; + +Physical Curve("ss2") = {5,6}; +Physical Curve("ss3") = {7,8}; +Physical Curve("bottom") = {2,3,4}; +Physical Curve("top") = {10,11}; +Physical Curve("right") = {9}; +Physical Curve("left") = {1}; + + + + diff --git a/src/meshmodifiers/InterfaceFromSideset.C b/src/meshmodifiers/InterfaceFromSideset.C index 066b6f49..6bc91426 100644 --- a/src/meshmodifiers/InterfaceFromSideset.C +++ b/src/meshmodifiers/InterfaceFromSideset.C @@ -24,7 +24,7 @@ validParams() InputParameters params = validParams(); params.addRequiredParam>( "sidesets", "The names of sidesets from which to create the new interface(s). " - "They must be straight and NOT have branches"); + "They MUST be straight and NOT have branches"); params.addRequiredParam>( "boundaries", "The names of sidesets forming the outside '" "boundaries of the whole mesh"); @@ -118,15 +118,38 @@ InterfaceFromSideset::modify() Elem * elem = (*it)->_elem; auto s = (*it)->_side; std::vector nodes_on_side = elem->nodes_on_side(s); - for (int n = 0; n < nodes_on_side.size(); ++n) + for (unsigned int n = 0; n < nodes_on_side.size(); ++n) boundary_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; + for (auto & sideset_name : sideset_names) + { + auto sidedset_id = _mesh_ptr->getBoundaryID(sideset_name); + 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) + { + auto current_node_id = elem->node_id(nodes_on_side[n]); + if (nb_neighbors_on_sideset[sideset_name].count(current_node_id) == 0) + nb_neighbors_on_sideset[sideset_name][current_node_id] = 1; + else + nb_neighbors_on_sideset[sideset_name][current_node_id] += 1; + } + } + } + // Find normal vector to plane of mesh by looking at first element auto normal_vec = getMeshNormalVector(mesh); // initialize the node to elemen map const auto & node_to_elem_map = _mesh_ptr->nodeToElemMap(); + std::map> new_node_to_elem_map; // loop on each provided sideset for (auto & sideset_name : sideset_names) @@ -136,24 +159,16 @@ InterfaceFromSideset::modify() std::set treated_node_ids; _new_boundary_sides_map.clear(); - // identify all nodes and nodes "not at the edge" of sideset + // identify all nodes on this sideset std::set sideset_node_ids; - std::map nb_neighbors_on_sideset; 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 (int n = 0; n < nodes_on_side.size(); ++n) - { - auto current_node_id = elem->node_id(nodes_on_side[n]); - sideset_node_ids.insert(current_node_id); - if (nb_neighbors_on_sideset.count(current_node_id) == 0) - nb_neighbors_on_sideset[current_node_id] = 1; - else - nb_neighbors_on_sideset[current_node_id] += 1; - } + for (unsigned int n = 0; n < nodes_on_side.size(); ++n) + sideset_node_ids.insert(elem->node_id(nodes_on_side[n])); } // loop on all elements involved in that sideset @@ -163,7 +178,7 @@ InterfaceFromSideset::modify() Elem * elem = (*it)->_elem; auto s = (*it)->_side; std::vector nodes_on_side = elem->nodes_on_side(s); - for (int n = 0; n < nodes_on_side.size(); ++n) + for (unsigned int n = 0; n < nodes_on_side.size(); ++n) { auto current_node_id = elem->node_id(nodes_on_side[n]); if (treated_node_ids.find(current_node_id) != treated_node_ids.end()) @@ -172,8 +187,8 @@ InterfaceFromSideset::modify() // Find out if we need to duplicate this node or not 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[current_node_id]>2; - //printf(" Node %d has %d neighbors on sideset\n", current_node_id, nb_neighbors_on_sideset[current_node_id]); + 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]); bool do_duplicate = false; if (is_node_on_boundary) do_duplicate = true; @@ -196,6 +211,11 @@ InterfaceFromSideset::modify() 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()); + new_node_to_elem_map[new_node->id()] = (node_to_elem_map.find(current_node_id))->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]; // Add sideset/boundary info to the new node std::vector node_boundary_ids = @@ -207,7 +227,11 @@ InterfaceFromSideset::modify() auto other_node_id = elem->node_id(nodes_on_side[(n+1)%2]); // retrieve connected elements from the map - auto node_to_elem_pair = node_to_elem_map.find(current_node_id); + std::map>::const_iterator node_to_elem_pair; + if (node_to_elem_map.find(current_node_id)==node_to_elem_map.end()) + node_to_elem_pair = new_node_to_elem_map.find(current_node_id); + else + node_to_elem_pair = node_to_elem_map.find(current_node_id); const std::vector & connected_elems = node_to_elem_pair->second; for (auto elem_id : connected_elems) @@ -239,18 +263,12 @@ InterfaceFromSideset::modify() { Elem * current_elem = mesh.elem_ptr(elem_id); Elem * connected_elem = mesh.elem_ptr(connected_elem_id); - //printf(" Testing pair %d - %d\n",elem_id, connected_elem_id); - //bool is_elem_on_that_side = elem_ids_on_that_side.find(elem_id) != elem_ids_on_that_side.end(); - //bool is_con_elem_on_other_side = elem_ids_on_other_side.find(connected_elem_id) != elem_ids_on_other_side.end(); - //printf(" is_elem_on_that_side=%d , is_con_elem_on_other_side=%d\n",is_elem_on_that_side, is_con_elem_on_other_side); if (current_elem != connected_elem && (elem_ids_on_that_side.find(elem_id) != elem_ids_on_that_side.end()) && (elem_ids_on_other_side.find(connected_elem_id) != elem_ids_on_other_side.end())) { - //printf(" selected pair\n"); if (current_elem->has_neighbor(connected_elem)) { - //printf(" add pair to mapping\n"); std::pair blocks_pair = std::make_pair(current_elem->subdomain_id(), connected_elem->subdomain_id()); _new_boundary_sides_map[blocks_pair].insert(std::make_pair( @@ -302,7 +320,7 @@ InterfaceFromSideset::isElementOnThatSideOfSegment( 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) { - int node_id = elem->node_id(i_v); + dof_id_type node_id = elem->node_id(i_v); bool is_node_on_sideset = (sideset_node_ids.find(node_id) != sideset_node_ids.end()); if ((node_id != node_id1) && (node_id != node_id2) && !is_node_on_sideset) { diff --git a/tests/meshmodifiers/gold/goldBreakMeshByBlock.i b/tests/meshmodifiers/gold/goldBreakMeshByBlock.i new file mode 100644 index 00000000..61a8e108 --- /dev/null +++ b/tests/meshmodifiers/gold/goldBreakMeshByBlock.i @@ -0,0 +1,9 @@ +[Mesh] + file = ../../meshes/fractureXtruncated.msh +[] + +[MeshModifiers] + [New_0] + type = BreakMeshByBlock + [] +[] diff --git a/tests/meshmodifiers/gold/goldBreakMeshByBlock_in.e b/tests/meshmodifiers/gold/goldBreakMeshByBlock_in.e new file mode 100644 index 00000000..9413ef8a Binary files /dev/null and b/tests/meshmodifiers/gold/goldBreakMeshByBlock_in.e differ diff --git a/tests/meshmodifiers/gold/testInterfaceFromSidesetDiffusion1.e b/tests/meshmodifiers/gold/testInterfaceFromSidesetDiffusion1.e new file mode 100644 index 00000000..c5868412 Binary files /dev/null and b/tests/meshmodifiers/gold/testInterfaceFromSidesetDiffusion1.e differ diff --git a/tests/meshmodifiers/gold/testInterfaceFromSidesetDiffusion2.e b/tests/meshmodifiers/gold/testInterfaceFromSidesetDiffusion2.e new file mode 100644 index 00000000..f583afc9 Binary files /dev/null and b/tests/meshmodifiers/gold/testInterfaceFromSidesetDiffusion2.e differ diff --git a/tests/meshmodifiers/gold/testInterfaceFromSidesetDiffusion3.e b/tests/meshmodifiers/gold/testInterfaceFromSidesetDiffusion3.e new file mode 100644 index 00000000..53a8b657 Binary files /dev/null and b/tests/meshmodifiers/gold/testInterfaceFromSidesetDiffusion3.e differ diff --git a/tests/meshmodifiers/interface_test.i b/tests/meshmodifiers/interface_test.i index e3e86818..72ee4084 100644 --- a/tests/meshmodifiers/interface_test.i +++ b/tests/meshmodifiers/interface_test.i @@ -1,6 +1,6 @@ [Mesh] type = FileMesh - file = test.msh + file = test2.msh [] [MeshModifiers] @@ -13,7 +13,7 @@ [break] type = InterfaceFromSideset boundaries = 'top bottom left right' - sidesets = 'ss3' + sidesets = 'ss2 ss3' split_interface = false [] [] @@ -36,14 +36,22 @@ [] [InterfaceKernels] - [interface] + [interface_ss2] type = InterfaceDarcy variable = p neighbor_var = p - boundary = interface + boundary = interface_ss2 fault_lewis_number = lewis_fault fault_thickness = 0.1 [] + [interface_ss3] + type = InterfaceDarcy + variable = p + neighbor_var = p + boundary = interface_ss3 + fault_lewis_number = lewis_fault + fault_thickness = 0.05 + [] [] [Materials] @@ -70,13 +78,13 @@ [left] type = DirichletBC variable = p - boundary = 'left_line' + boundary = 'left' value = 1 [] [right] type = DirichletBC variable = p - boundary = 'right_line' + boundary = 'right' value = 0 [] [] diff --git a/tests/meshmodifiers/testInterfaceFromSideset.i b/tests/meshmodifiers/testInterfaceFromSideset.i new file mode 100644 index 00000000..8680915c --- /dev/null +++ b/tests/meshmodifiers/testInterfaceFromSideset.i @@ -0,0 +1,14 @@ +# Test for InterfaceFromSideset MeshModifier on 2D block +# with simple X fractures (touching all sides) + +[Mesh] + file = ../../meshes/fractureX.msh +[] + +[MeshModifiers] + [New_0] + type = InterfaceFromSideset + sidesets = 'ss2 ss3' + boundaries = 'left right top bottom' + [] +[] diff --git a/tests/meshmodifiers/testInterfaceFromSidesetDiffusion1.i b/tests/meshmodifiers/testInterfaceFromSidesetDiffusion1.i new file mode 100644 index 00000000..f623fe61 --- /dev/null +++ b/tests/meshmodifiers/testInterfaceFromSidesetDiffusion1.i @@ -0,0 +1,50 @@ +[Mesh] + type = FileMesh + file = ../../meshes/fractureX.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 = testInterfaceFromSidesetDiffusion1 + exodus = true +[] diff --git a/tests/meshmodifiers/testInterfaceFromSidesetDiffusion2.i b/tests/meshmodifiers/testInterfaceFromSidesetDiffusion2.i new file mode 100644 index 00000000..bc6fc9d5 --- /dev/null +++ b/tests/meshmodifiers/testInterfaceFromSidesetDiffusion2.i @@ -0,0 +1,50 @@ +[Mesh] + type = FileMesh + file = ../../meshes/fractureX.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 = testInterfaceFromSidesetDiffusion2 + exodus = true +[] diff --git a/tests/meshmodifiers/testInterfaceFromSidesetDiffusion3.i b/tests/meshmodifiers/testInterfaceFromSidesetDiffusion3.i new file mode 100644 index 00000000..528000c6 --- /dev/null +++ b/tests/meshmodifiers/testInterfaceFromSidesetDiffusion3.i @@ -0,0 +1,50 @@ +[Mesh] + type = FileMesh + file = ../../meshes/fractureXincluded.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 = testInterfaceFromSidesetDiffusion3 + exodus = true +[] diff --git a/tests/meshmodifiers/tests b/tests/meshmodifiers/tests new file mode 100644 index 00000000..db9d9394 --- /dev/null +++ b/tests/meshmodifiers/tests @@ -0,0 +1,27 @@ +[Tests] + [./test_interfaceFromSideset] + # Should produce the same mesh as BreakMeshByBlock meshmodifier + type = 'Exodiff' + input = 'testinterfaceFromSideset.i' + exodiff = 'goldBreakMeshByBlock_in.e' + cli_args = '--mesh-only=goldBreakMeshByBlock_in.e' + [../] + [./test_interfaceFromSideset_diff1] + # Testing diffusion on X reaching sides from left to right + type = 'Exodiff' + input = 'testInterfaceFromSidesetDiffusion1.i' + exodiff = 'testInterfaceFromSidesetDiffusion1.e' + [../] + [./test_interfaceFromSideset_diff2] + # Testing diffusion on X reaching sides from right to left + type = 'Exodiff' + input = 'testInterfaceFromSidesetDiffusion2.i' + exodiff = 'testInterfaceFromSidesetDiffusion2.e' + [../] + [./test_interfaceFromSideset_diff3] + # Testing diffusion on X not reaching sides + type = 'Exodiff' + input = 'testInterfaceFromSidesetDiffusion3.i' + exodiff = 'testInterfaceFromSidesetDiffusion3.e' + [../] +[]