Skip to content

Commit

Permalink
#180 supporting T junctions
Browse files Browse the repository at this point in the history
  • Loading branch information
@pou036
pou036 committed Aug 9, 2019
1 parent a5267d8 commit 41792fd
Show file tree
Hide file tree
Showing 9 changed files with 321 additions and 27 deletions.
67 changes: 41 additions & 26 deletions src/meshmodifiers/InterfaceFromSideset.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<int> boundary_node_ids;
std::set<int> allsidesets_node_ids;
std::map<BoundaryName, std::set<int>> 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;
Expand All @@ -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<unsigned int> 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<BoundaryName, std::map<int,int>> nb_neighbors_on_sideset;
Expand Down Expand Up @@ -118,17 +137,12 @@ InterfaceFromSideset::modify()
std::set<int> treated_node_ids;
_new_boundary_sides_map.clear();

// identify all nodes on this sideset
std::set<int> 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<unsigned int> 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<int> 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)
Expand All @@ -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<BoundaryID> elem_ids_on_that_side;
std::set<BoundaryID> elem_ids_on_other_side;
Expand All @@ -166,19 +183,26 @@ 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<dof_id_type, std::vector<dof_id_type>>::const_iterator
node_to_elem_pair = node_to_elem_map.find(current_node_id);
if (node_to_elem_pair==node_to_elem_map.end())
node_to_elem_pair = new_node_to_elem_map.find(current_node_id);
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<boundary_id_type> node_boundary_ids =
Expand All @@ -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);
Expand Down Expand Up @@ -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)
{
Expand All @@ -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");
Expand Down Expand Up @@ -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;
}

Expand All @@ -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;
}

Expand Down
2 changes: 1 addition & 1 deletion tests/meshmodifiers/fractureX.geo
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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};

Expand Down
52 changes: 52 additions & 0 deletions tests/meshmodifiers/fractureY.geo
View file
Original file line number Diff line number Diff line change
@@ -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};




115 changes: 115 additions & 0 deletions tests/meshmodifiers/fractureY.msh
View file
Original file line number Diff line number Diff line change
@@ -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
Binary file added tests/meshmodifiers/gold/testInterfaceFromSidesetDiffusion4.e
View file
Binary file not shown.
Binary file added tests/meshmodifiers/gold/testInterfaceFromSidesetDiffusion5.e
View file
Binary file not shown.
50 changes: 50 additions & 0 deletions tests/meshmodifiers/testInterfaceFromSidesetDiffusion4.i
View file
Original file line number Diff line number Diff line change
@@ -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
[]
Loading

0 comments on commit 41792fd

Please sign in to comment.