Skip to content

Commit

Permalink
#180 first version running, along with tests (but bad_alloc problem a…
Browse files Browse the repository at this point in the history
…ppearing sometimes, still needs debugging...)
  • Loading branch information
pou036 committed Aug 9, 2019
1 parent 542dddd commit 3d34974
Show file tree
Hide file tree
Showing 16 changed files with 403 additions and 31 deletions.
1 change: 0 additions & 1 deletion include/meshmodifiers/InterfaceFromSideset.h
Original file line number Diff line number Diff line change
Expand Up @@ -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<int>, const std::vector<Real>);
std::vector<Real> getMeshNormalVector(const MeshBase &);

std::set<std::pair<subdomain_id_type, subdomain_id_type>> _neighboring_block_list;
std::map<std::pair<subdomain_id_type, subdomain_id_type>,
std::set<std::pair<dof_id_type, unsigned int>>>
_new_boundary_sides_map;
Expand Down
58 changes: 58 additions & 0 deletions meshes/fractureX.geo
Original file line number Diff line number Diff line change
@@ -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};




43 changes: 43 additions & 0 deletions meshes/fractureXincluded.geo
Original file line number Diff line number Diff line change
@@ -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};




46 changes: 46 additions & 0 deletions meshes/fractureXtruncated.geo
Original file line number Diff line number Diff line change
@@ -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};




66 changes: 42 additions & 24 deletions src/meshmodifiers/InterfaceFromSideset.C
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ validParams<InterfaceFromSideset>()
InputParameters params = validParams<BreakMeshByBlockBase>();
params.addRequiredParam<std::vector<BoundaryName>>(
"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<std::vector<BoundaryName>>(
"boundaries", "The names of sidesets forming the outside '"
"boundaries of the whole mesh");
Expand Down Expand Up @@ -118,15 +118,38 @@ InterfaceFromSideset::modify()
Elem * elem = (*it)->_elem;
auto s = (*it)->_side;
std::vector<unsigned int> 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<BoundaryName, std::map<int,int>> 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<unsigned int> 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<dof_id_type, std::vector<dof_id_type>> new_node_to_elem_map;

// loop on each provided sideset
for (auto & sideset_name : sideset_names)
Expand All @@ -136,24 +159,16 @@ InterfaceFromSideset::modify()
std::set<int> 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<int> sideset_node_ids;
std::map<int,int> 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<unsigned int> 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
Expand All @@ -163,7 +178,7 @@ InterfaceFromSideset::modify()
Elem * elem = (*it)->_elem;
auto s = (*it)->_side;
std::vector<unsigned int> 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())
Expand All @@ -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;
Expand All @@ -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<boundary_id_type> node_boundary_ids =
Expand All @@ -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<dof_id_type, std::vector<dof_id_type>>::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<dof_id_type> & connected_elems =
node_to_elem_pair->second;
for (auto elem_id : connected_elems)
Expand Down Expand Up @@ -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<subdomain_id_type, subdomain_id_type> blocks_pair =
std::make_pair(current_elem->subdomain_id(), connected_elem->subdomain_id());
_new_boundary_sides_map[blocks_pair].insert(std::make_pair(
Expand Down Expand Up @@ -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)
{
Expand Down
9 changes: 9 additions & 0 deletions tests/meshmodifiers/gold/goldBreakMeshByBlock.i
Original file line number Diff line number Diff line change
@@ -0,0 +1,9 @@
[Mesh]
file = ../../meshes/fractureXtruncated.msh
[]

[MeshModifiers]
[New_0]
type = BreakMeshByBlock
[]
[]
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
20 changes: 14 additions & 6 deletions tests/meshmodifiers/interface_test.i
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
[Mesh]
type = FileMesh
file = test.msh
file = test2.msh
[]

[MeshModifiers]
Expand All @@ -13,7 +13,7 @@
[break]
type = InterfaceFromSideset
boundaries = 'top bottom left right'
sidesets = 'ss3'
sidesets = 'ss2 ss3'
split_interface = false
[]
[]
Expand All @@ -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]
Expand All @@ -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
[]
[]
Expand Down
14 changes: 14 additions & 0 deletions tests/meshmodifiers/testInterfaceFromSideset.i
Original file line number Diff line number Diff line change
@@ -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'
[]
[]
Loading

0 comments on commit 3d34974

Please sign in to comment.