Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F80276634
cohesive_element_inserter_helper.cc
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Fri, Aug 30, 21:28
Size
32 KB
Mime Type
text/x-c
Expires
Sun, Sep 1, 21:28 (1 d, 22 h)
Engine
blob
Format
Raw Data
Handle
20332397
Attached To
rAKA akantu
cohesive_element_inserter_helper.cc
View Options
/**
* @file cohesive_element_inserter_helper.cc
*
* @author Nicolas Richart
*
* @date creation jeu sep 03 2020
*
* @brief A Documented file.
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "cohesive_element_inserter_helper.hh"
/* -------------------------------------------------------------------------- */
#include "data_accessor.hh"
#include "element_synchronizer.hh"
#include "fe_engine.hh"
#include "mesh_accessor.hh"
/* -------------------------------------------------------------------------- */
#include <queue>
/* -------------------------------------------------------------------------- */
namespace
akantu
{
/* -------------------------------------------------------------------------- */
CohesiveElementInserterHelper
::
CohesiveElementInserterHelper
(
Mesh
&
mesh
,
const
ElementTypeMapArray
<
bool
>
&
facet_insertion
)
:
doubled_nodes
(
0
,
2
,
"doubled_nodes"
),
mesh
(
mesh
),
mesh_facets
(
mesh
.
getMeshFacets
())
{
auto
spatial_dimension
=
mesh_facets
.
getSpatialDimension
();
for
(
auto
gt
:
ghost_types
)
{
for
(
auto
type
:
mesh_facets
.
elementTypes
(
_ghost_type
=
gt
))
{
nb_new_facets
(
type
,
gt
)
=
mesh_facets
.
getNbElement
(
type
,
gt
);
}
}
std
::
array
<
Int
,
2
>
nb_facet_to_insert
{
0
,
0
};
// creates a vector of the facets to insert
std
::
vector
<
Element
>
potential_facets_to_double
;
for
(
auto
gt_facet
:
ghost_types
)
{
for
(
auto
type_facet
:
mesh_facets
.
elementTypes
(
spatial_dimension
-
1
,
gt_facet
))
{
const
auto
&
f_insertion
=
facet_insertion
(
type_facet
,
gt_facet
);
auto
&
counter
=
nb_facet_to_insert
[
gt_facet
==
_not_ghost
?
0
:
1
];
for
(
auto
&&
data
:
enumerate
(
f_insertion
))
{
if
(
std
::
get
<
1
>
(
data
))
{
UInt
el
=
std
::
get
<
0
>
(
data
);
potential_facets_to_double
.
push_back
({
type_facet
,
el
,
gt_facet
});
++
counter
;
}
}
}
}
// Defines a global order of insertion oof cohesive element to ensure node
// doubling appens in the smae order, this is necessary for the global node
// numbering
if
(
mesh_facets
.
isDistributed
())
{
const
auto
&
comm
=
mesh_facets
.
getCommunicator
();
ElementTypeMapArray
<
Int
>
global_orderings
;
global_orderings
.
initialize
(
mesh_facets
,
_spatial_dimension
=
spatial_dimension
-
1
,
_with_nb_element
=
true
,
_default_value
=
-
1
);
auto
starting_index
=
nb_facet_to_insert
[
0
];
comm
.
exclusiveScan
(
starting_index
);
// define the local numbers for facet to insert
for
(
auto
gt_facet
:
ghost_types
)
{
for
(
auto
type_facet
:
mesh_facets
.
elementTypes
(
spatial_dimension
-
1
,
gt_facet
))
{
for
(
auto
data
:
zip
(
facet_insertion
(
type_facet
,
gt_facet
),
global_orderings
(
type_facet
,
gt_facet
)))
{
if
(
std
::
get
<
0
>
(
data
))
{
std
::
get
<
1
>
(
data
)
=
starting_index
;
++
starting_index
;
}
}
}
}
// retreives the oorder number from neighoring processors
auto
&&
synchronizer
=
mesh_facets
.
getElementSynchronizer
();
SimpleElementDataAccessor
<
Int
>
data_accessor
(
global_orderings
,
SynchronizationTag
::
_ce_insertion_order
);
synchronizer
.
synchronizeOnce
(
data_accessor
,
SynchronizationTag
::
_ce_insertion_order
);
// sort the facets to double based on this global ordering
std
::
sort
(
potential_facets_to_double
.
begin
(),
potential_facets_to_double
.
end
(),
[
&
global_orderings
](
auto
&&
el1
,
auto
&&
el2
)
{
return
global_orderings
(
el1
)
<
global_orderings
(
el2
);
});
}
for
(
auto
d
:
arange
(
spatial_dimension
))
{
facets_to_double_by_dim
[
d
]
=
std
::
make_unique
<
Array
<
Element
>>
(
0
,
2
,
"facets_to_double_"
+
std
::
to_string
(
d
));
}
auto
&
facets_to_double
=
*
facets_to_double_by_dim
[
spatial_dimension
-
1
];
MeshAccessor
mesh_accessor
(
mesh_facets
);
auto
&
elements_to_subelements
=
mesh_accessor
.
getElementToSubelement
();
for
(
auto
&&
facet_to_double
:
potential_facets_to_double
)
{
auto
gt_facet
=
facet_to_double
.
ghost_type
;
auto
type_facet
=
facet_to_double
.
type
;
auto
&
elements_to_facets
=
elements_to_subelements
(
type_facet
,
gt_facet
);
auto
&
elements_to_facet
=
elements_to_facets
(
facet_to_double
.
element
);
if
(
elements_to_facet
[
1
].
type
==
_not_defined
#if defined(AKANTU_COHESIVE_ELEMENT)
||
elements_to_facet
[
1
].
kind
()
==
_ek_cohesive
#endif
)
{
AKANTU_DEBUG_WARNING
(
"attempt to double a facet on the boundary"
);
continue
;
}
auto
new_facet
=
nb_new_facets
(
type_facet
,
gt_facet
)
++
;
facets_to_double
.
push_back
(
Vector
<
Element
>
{
facet_to_double
,
Element
{
type_facet
,
new_facet
,
gt_facet
}});
/// update facet_to_element vector
auto
&
element_to_update
=
elements_to_facet
[
1
];
auto
el
=
element_to_update
.
element
;
auto
&
facets_to_elements
=
mesh_facets
.
getSubelementToElement
(
element_to_update
.
type
,
element_to_update
.
ghost_type
);
auto
facets_to_element
=
Vector
<
Element
>
(
make_view
(
facets_to_elements
,
facets_to_elements
.
getNbComponent
())
.
begin
()[
el
]);
auto
facet_to_update
=
std
::
find
(
facets_to_element
.
begin
(),
facets_to_element
.
end
(),
facet_to_double
);
AKANTU_DEBUG_ASSERT
(
facet_to_update
!=
facets_to_element
.
end
(),
"Facet not found"
);
facet_to_update
->
element
=
new_facet
;
/// update elements connected to facet
const
auto
&
first_facet_list
=
elements_to_facet
;
elements_to_facets
.
push_back
(
first_facet_list
);
/// set new and original facets as boundary facets
elements_to_facets
(
new_facet
)[
0
]
=
elements_to_facets
(
new_facet
)[
1
];
elements_to_facets
(
new_facet
)[
1
]
=
ElementNull
;
elements_to_facets
(
facet_to_double
.
element
)[
1
]
=
ElementNull
;
}
}
/* -------------------------------------------------------------------------- */
inline
auto
CohesiveElementInserterHelper
::
hasElement
(
const
Vector
<
UInt
>
&
nodes_element
,
const
Vector
<
UInt
>
&
nodes
)
->
bool
{
// if one of the nodes of "nodes" is not in "nodes_element" it stops
auto
it
=
std
::
mismatch
(
nodes
.
begin
(),
nodes
.
end
(),
nodes_element
.
begin
(),
[
&
](
auto
&&
node
,
auto
&&
/*node2*/
)
->
bool
{
auto
it
=
std
::
find
(
nodes_element
.
begin
(),
nodes_element
.
end
(),
node
);
return
(
it
!=
nodes_element
.
end
());
});
// true if all "nodes" where found in "nodes_element"
return
(
it
.
first
==
nodes
.
end
());
}
/* -------------------------------------------------------------------------- */
inline
auto
CohesiveElementInserterHelper
::
removeElementsInVector
(
const
std
::
vector
<
Element
>
&
elem_to_remove
,
std
::
vector
<
Element
>
&
elem_list
)
->
bool
{
if
(
elem_list
.
size
()
<=
elem_to_remove
.
size
())
{
return
true
;
}
auto
el_it
=
elem_to_remove
.
begin
();
auto
el_last
=
elem_to_remove
.
end
();
std
::
vector
<
Element
>::
iterator
el_del
;
UInt
deletions
=
0
;
for
(;
el_it
!=
el_last
;
++
el_it
)
{
el_del
=
std
::
find
(
elem_list
.
begin
(),
elem_list
.
end
(),
*
el_it
);
if
(
el_del
!=
elem_list
.
end
())
{
elem_list
.
erase
(
el_del
);
++
deletions
;
}
}
AKANTU_DEBUG_ASSERT
(
deletions
==
0
||
deletions
==
elem_to_remove
.
size
(),
"Not all elements have been erased"
);
return
deletions
==
0
;
}
/* -------------------------------------------------------------------------- */
void
CohesiveElementInserterHelper
::
updateElementalConnectivity
(
Mesh
&
mesh
,
UInt
old_node
,
UInt
new_node
,
const
std
::
vector
<
Element
>
&
element_list
,
const
std
::
vector
<
Element
>
*
facet_list
)
{
AKANTU_DEBUG_IN
();
for
(
auto
&
element
:
element_list
)
{
if
(
element
.
type
==
_not_defined
)
{
continue
;
}
Vector
<
UInt
>
connectivity
=
mesh
.
getConnectivity
(
element
);
if
(
element
.
kind
()
==
_ek_cohesive
)
{
AKANTU_DEBUG_ASSERT
(
facet_list
!=
nullptr
,
"Provide a facet list in order to update cohesive elements"
);
const
Vector
<
Element
>
facets
=
mesh_facets
.
getSubelementToElement
(
element
);
auto
facet_nb_nodes
=
connectivity
.
size
()
/
2
;
/// loop over cohesive element's facets
for
(
const
auto
&
facet
:
enumerate
(
facets
))
{
/// skip facets if not present in the list
if
(
std
::
find
(
facet_list
->
begin
(),
facet_list
->
end
(),
std
::
get
<
1
>
(
facet
))
==
facet_list
->
end
())
{
continue
;
}
auto
n
=
std
::
get
<
0
>
(
facet
);
auto
begin
=
connectivity
.
begin
()
+
static_cast
<
UInt
>
(
n
*
facet_nb_nodes
);
auto
end
=
begin
+
facet_nb_nodes
;
auto
it
=
std
::
find
(
begin
,
end
,
old_node
);
AKANTU_DEBUG_ASSERT
(
it
!=
end
,
"Node not found in current element"
);
*
it
=
new_node
;
}
}
else
{
auto
it
=
std
::
find
(
connectivity
.
begin
(),
connectivity
.
end
(),
old_node
);
AKANTU_DEBUG_ASSERT
(
it
!=
connectivity
.
end
(),
"Node not found in current element"
);
/// update connectivity
*
it
=
new_node
;
}
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
CohesiveElementInserterHelper
::
updateSubelementToElement
(
UInt
dim
,
bool
facet_mode
)
{
auto
&
facets_to_double
=
*
facets_to_double_by_dim
[
dim
];
auto
&
facets_to_subfacets
=
elementsOfDimToElementsOfDim
(
dim
+
static_cast
<
decltype
(
dim
)
>
(
facet_mode
),
dim
);
for
(
auto
&&
data
:
zip
(
make_view
(
facets_to_double
,
2
),
facets_to_subfacets
))
{
const
auto
&
old_subfacet
=
std
::
get
<
0
>
(
data
)[
0
];
const
auto
&
new_subfacet
=
std
::
get
<
0
>
(
data
)[
1
];
auto
&
facet_to_subfacets
=
std
::
get
<
1
>
(
data
);
MeshAccessor
mesh_accessor
(
mesh_facets
);
for
(
auto
&
facet
:
facet_to_subfacets
)
{
Vector
<
Element
>
subfacets
=
mesh_accessor
.
getSubelementToElement
(
facet
);
auto
&&
subfacet_to_update_it
=
std
::
find
(
subfacets
.
begin
(),
subfacets
.
end
(),
old_subfacet
);
AKANTU_DEBUG_ASSERT
(
subfacet_to_update_it
!=
subfacets
.
end
(),
"Subfacet not found"
);
*
subfacet_to_update_it
=
new_subfacet
;
}
}
}
/* -------------------------------------------------------------------------- */
void
CohesiveElementInserterHelper
::
updateElementToSubelement
(
UInt
dim
,
bool
facet_mode
)
{
auto
&
facets_to_double
=
*
facets_to_double_by_dim
[
dim
];
auto
&
facets_to_subfacets
=
elementsOfDimToElementsOfDim
(
dim
+
static_cast
<
decltype
(
dim
)
>
(
facet_mode
),
dim
);
MeshAccessor
mesh_accessor
(
mesh_facets
);
// resize the arrays
mesh_accessor
.
getElementToSubelement
().
initialize
(
mesh_facets
,
_spatial_dimension
=
dim
,
_with_nb_element
=
true
);
for
(
auto
&&
data
:
zip
(
make_view
(
facets_to_double
,
2
),
facets_to_subfacets
))
{
const
auto
&
new_facet
=
std
::
get
<
0
>
(
data
)[
1
];
mesh_accessor
.
getElementToSubelement
(
new_facet
)
=
std
::
get
<
1
>
(
data
);
}
}
/* -------------------------------------------------------------------------- */
void
CohesiveElementInserterHelper
::
updateQuadraticSegments
(
UInt
dim
)
{
AKANTU_DEBUG_IN
();
auto
spatial_dimension
=
mesh
.
getSpatialDimension
();
auto
&
facets_to_double
=
*
facets_to_double_by_dim
[
dim
];
MeshAccessor
mesh_accessor
(
mesh_facets
);
auto
&
connectivities
=
mesh_accessor
.
getConnectivities
();
/// this ones matter only for segments in 3D
Array
<
std
::
vector
<
Element
>>
*
element_to_subfacet_double
=
nullptr
;
Array
<
std
::
vector
<
Element
>>
*
facet_to_subfacet_double
=
nullptr
;
if
(
dim
==
spatial_dimension
-
2
)
{
element_to_subfacet_double
=
&
elementsOfDimToElementsOfDim
(
dim
+
2
,
dim
);
facet_to_subfacet_double
=
&
elementsOfDimToElementsOfDim
(
dim
+
1
,
dim
);
}
auto
&
element_to_subelement
=
mesh_facets
.
getElementToSubelement
();
std
::
vector
<
UInt
>
middle_nodes
;
for
(
auto
&&
facet_to_double
:
make_view
(
facets_to_double
,
2
))
{
auto
&
old_facet
=
facet_to_double
[
0
];
if
(
old_facet
.
type
!=
_segment_3
)
{
continue
;
}
auto
node
=
connectivities
(
old_facet
,
2
);
if
(
not
mesh
.
isPureGhostNode
(
node
))
{
middle_nodes
.
push_back
(
node
);
}
}
auto
n
=
doubled_nodes
.
size
();
doubleNodes
(
middle_nodes
);
for
(
auto
&&
data
:
enumerate
(
make_view
(
facets_to_double
,
2
)))
{
auto
facet
=
std
::
get
<
0
>
(
data
);
auto
&
old_facet
=
std
::
get
<
1
>
(
data
)[
0
];
if
(
old_facet
.
type
!=
_segment_3
)
{
continue
;
}
auto
old_node
=
connectivities
(
old_facet
,
2
);
if
(
mesh
.
isPureGhostNode
(
old_node
))
{
continue
;
}
auto
new_node
=
doubled_nodes
(
n
,
1
);
auto
&
new_facet
=
std
::
get
<
1
>
(
data
)[
1
];
connectivities
(
new_facet
,
2
)
=
new_node
;
if
(
dim
==
spatial_dimension
-
2
)
{
updateElementalConnectivity
(
mesh_facets
,
old_node
,
new_node
,
element_to_subelement
(
new_facet
,
0
));
updateElementalConnectivity
(
mesh
,
old_node
,
new_node
,
(
*
element_to_subfacet_double
)(
facet
),
&
(
*
facet_to_subfacet_double
)(
facet
));
}
else
{
updateElementalConnectivity
(
mesh
,
old_node
,
new_node
,
element_to_subelement
(
new_facet
,
0
));
}
++
n
;
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
UInt
CohesiveElementInserterHelper
::
insertCohesiveElement
()
{
auto
nb_new_facets
=
insertFacetsOnly
();
if
(
nb_new_facets
==
0
)
{
return
0
;
}
updateCohesiveData
();
return
nb_new_facets
;
}
/* -------------------------------------------------------------------------- */
UInt
CohesiveElementInserterHelper
::
insertFacetsOnly
()
{
UInt
spatial_dimension
=
mesh
.
getSpatialDimension
();
if
(
facets_to_double_by_dim
[
spatial_dimension
-
1
]
->
empty
())
{
return
0
;
}
if
(
spatial_dimension
==
1
)
{
doublePointFacet
();
}
else
if
(
spatial_dimension
==
2
)
{
doubleFacets
<
1
>
();
findSubfacetToDouble
<
1
>
();
doubleSubfacet
<
2
>
();
}
else
if
(
spatial_dimension
==
3
)
{
doubleFacets
<
2
>
();
findSubfacetToDouble
<
2
>
();
doubleFacets
<
1
>
();
findSubfacetToDouble
<
1
>
();
doubleSubfacet
<
3
>
();
}
return
facets_to_double_by_dim
[
spatial_dimension
-
1
]
->
size
();
}
/* -------------------------------------------------------------------------- */
template
<
UInt
dim
>
void
CohesiveElementInserterHelper
::
doubleFacets
()
{
AKANTU_DEBUG_IN
();
NewElementsEvent
new_facets
;
auto
spatial_dimension
=
mesh_facets
.
getSpatialDimension
();
auto
&
facets_to_double
=
*
facets_to_double_by_dim
[
dim
];
MeshAccessor
mesh_accessor
(
mesh_facets
);
for
(
auto
&&
facet_to_double
:
make_view
(
facets_to_double
,
2
))
{
auto
&&
old_facet
=
facet_to_double
[
0
];
auto
&&
new_facet
=
facet_to_double
[
1
];
auto
&
facets_connectivities
=
mesh_accessor
.
getConnectivity
(
old_facet
.
type
,
old_facet
.
ghost_type
);
facets_connectivities
.
resize
(
nb_new_facets
(
old_facet
.
type
,
old_facet
.
ghost_type
));
auto
facets_connectivities_begin
=
make_view
(
facets_connectivities
,
facets_connectivities
.
getNbComponent
())
.
begin
();
// copy the connectivities
Vector
<
UInt
>
new_conn
(
facets_connectivities_begin
[
new_facet
.
element
]);
Vector
<
UInt
>
old_conn
(
facets_connectivities_begin
[
old_facet
.
element
]);
new_conn
=
old_conn
;
// this will fail if multiple facet types exists for a given element type
// \TODO handle multiple sub-facet types
auto
nb_subfacet_per_facet
=
Mesh
::
getNbFacetsPerElement
(
old_facet
.
type
);
auto
&
subfacets_to_facets
=
mesh_accessor
.
getSubelementToElementNC
(
old_facet
.
type
,
old_facet
.
ghost_type
);
subfacets_to_facets
.
resize
(
nb_new_facets
(
old_facet
.
type
,
old_facet
.
ghost_type
),
ElementNull
);
auto
subfacets_to_facets_begin
=
make_view
(
subfacets_to_facets
,
nb_subfacet_per_facet
).
begin
();
// copy the subfacet to facets information
Vector
<
Element
>
old_subfacets_to_facet
(
subfacets_to_facets_begin
[
old_facet
.
element
]);
Vector
<
Element
>
new_subfacet_to_facet
(
subfacets_to_facets_begin
[
new_facet
.
element
]);
new_subfacet_to_facet
=
old_subfacets_to_facet
;
for
(
auto
&
subfacet
:
old_subfacets_to_facet
)
{
if
(
subfacet
==
ElementNull
)
{
continue
;
}
/// update facet_to_subfacet array
mesh_accessor
.
getElementToSubelement
(
subfacet
).
push_back
(
new_facet
);
}
new_facets
.
getList
().
push_back
(
new_facet
);
}
/// update facet_to_subfacet and _segment_3 facets if any
if
(
dim
!=
spatial_dimension
-
1
)
{
updateSubelementToElement
(
dim
,
true
);
updateElementToSubelement
(
dim
,
true
);
updateQuadraticSegments
(
dim
);
}
else
{
updateQuadraticSegments
(
dim
);
}
mesh_facets
.
sendEvent
(
new_facets
);
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
UInt
dim
>
void
CohesiveElementInserterHelper
::
findSubfacetToDouble
()
{
AKANTU_DEBUG_IN
();
UInt
spatial_dimension
=
mesh_facets
.
getSpatialDimension
();
MeshAccessor
mesh_accessor
(
mesh_facets
);
auto
&
facets_to_double
=
*
facets_to_double_by_dim
[
spatial_dimension
-
1
];
auto
&
subfacets_to_facets
=
mesh_facets
.
getSubelementToElement
();
auto
&
elements_to_facets
=
mesh_accessor
.
getElementToSubelement
();
/// loop on every facet
for
(
auto
f
:
arange
(
2
))
{
for
(
auto
&&
facet_to_double
:
make_view
(
facets_to_double
,
2
))
{
auto
old_facet
=
facet_to_double
[
f
];
auto
new_facet
=
facet_to_double
[
1
-
f
];
const
auto
&
starting_element
=
elements_to_facets
(
new_facet
,
0
)[
0
];
auto
current_facet
=
old_facet
;
Vector
<
Element
>
subfacets_to_facet
=
subfacets_to_facets
.
get
(
old_facet
);
/// loop on every subfacet
for
(
auto
&
subfacet
:
subfacets_to_facet
)
{
if
(
subfacet
==
ElementNull
)
{
continue
;
}
if
(
dim
==
spatial_dimension
-
2
)
{
Vector
<
Element
>
subsubfacets_to_subfacet
=
subfacets_to_facets
.
get
(
subfacet
);
/// loop on every subsubfacet
for
(
auto
&
subsubfacet
:
subsubfacets_to_subfacet
)
{
if
(
subsubfacet
==
ElementNull
)
{
continue
;
}
Vector
<
UInt
>
subsubfacet_connectivity
=
mesh_facets
.
getConnectivity
(
subsubfacet
);
std
::
vector
<
Element
>
element_list
;
std
::
vector
<
Element
>
facet_list
;
std
::
vector
<
Element
>
subfacet_list
;
/// check if subsubfacet is to be doubled
if
(
findElementsAroundSubfacet
(
starting_element
,
current_facet
,
subsubfacet_connectivity
,
element_list
,
facet_list
,
&
subfacet_list
)
==
false
and
removeElementsInVector
(
subfacet_list
,
elements_to_facets
(
subsubfacet
))
==
false
)
{
Element
new_subsubfacet
{
subsubfacet
.
type
,
nb_new_facets
(
subsubfacet
.
type
,
subsubfacet
.
ghost_type
)
++
,
subsubfacet
.
ghost_type
};
facets_to_double_by_dim
[
dim
-
1
]
->
push_back
(
Vector
<
Element
>
{
subsubfacet
,
new_subsubfacet
});
elementsOfDimToElementsOfDim
(
dim
-
1
,
dim
-
1
)
.
push_back
(
subfacet_list
);
elementsOfDimToElementsOfDim
(
dim
,
dim
-
1
)
.
push_back
(
facet_list
);
elementsOfDimToElementsOfDim
(
dim
+
1
,
dim
-
1
)
.
push_back
(
element_list
);
}
}
}
else
{
std
::
vector
<
Element
>
element_list
;
std
::
vector
<
Element
>
facet_list
;
Vector
<
UInt
>
subfacet_connectivity
=
mesh_facets
.
getConnectivity
(
subfacet
);
/// check if subfacet is to be doubled
if
(
findElementsAroundSubfacet
(
starting_element
,
current_facet
,
subfacet_connectivity
,
element_list
,
facet_list
)
==
false
and
removeElementsInVector
(
facet_list
,
elements_to_facets
(
subfacet
))
==
false
)
{
Element
new_subfacet
{
subfacet
.
type
,
nb_new_facets
(
subfacet
.
type
,
subfacet
.
ghost_type
)
++
,
subfacet
.
ghost_type
};
facets_to_double_by_dim
[
dim
-
1
]
->
push_back
(
Vector
<
Element
>
{
subfacet
,
new_subfacet
});
elementsOfDimToElementsOfDim
(
dim
,
dim
-
1
).
push_back
(
facet_list
);
elementsOfDimToElementsOfDim
(
dim
+
1
,
dim
-
1
)
.
push_back
(
element_list
);
}
}
}
}
}
}
/* -------------------------------------------------------------------------- */
void
CohesiveElementInserterHelper
::
doubleNodes
(
const
std
::
vector
<
UInt
>
&
old_nodes
)
{
auto
&
position
=
mesh
.
getNodes
();
auto
spatial_dimension
=
mesh
.
getSpatialDimension
();
auto
old_nb_nodes
=
position
.
size
();
position
.
reserve
(
old_nb_nodes
+
old_nodes
.
size
());
doubled_nodes
.
reserve
(
doubled_nodes
.
size
()
+
old_nodes
.
size
());
auto
position_begin
=
position
.
begin
(
spatial_dimension
);
for
(
auto
&&
data
:
enumerate
(
old_nodes
))
{
auto
n
=
std
::
get
<
0
>
(
data
);
auto
old_node
=
std
::
get
<
1
>
(
data
);
decltype
(
old_node
)
new_node
=
old_nb_nodes
+
n
;
/// store doubled nodes
doubled_nodes
.
push_back
(
Vector
<
UInt
>
{
old_node
,
new_node
});
/// update position
Vector
<
Real
>
coords
=
Vector
<
Real
>
(
position_begin
[
old_node
]);
position
.
push_back
(
coords
);
}
}
/* -------------------------------------------------------------------------- */
bool
CohesiveElementInserterHelper
::
findElementsAroundSubfacet
(
const
Element
&
starting_element
,
const
Element
&
end_facet
,
const
Vector
<
UInt
>
&
subfacet_connectivity
,
std
::
vector
<
Element
>
&
element_list
,
std
::
vector
<
Element
>
&
facet_list
,
std
::
vector
<
Element
>
*
subfacet_list
)
{
bool
facet_matched
=
false
;
element_list
.
push_back
(
starting_element
);
std
::
queue
<
Element
>
elements_to_check
;
elements_to_check
.
push
(
starting_element
);
/// keep going as long as there are elements to check
while
(
not
elements_to_check
.
empty
())
{
/// check current element
Element
&
current_element
=
elements_to_check
.
front
();
const
Vector
<
Element
>
facets_to_element
=
mesh_facets
.
getSubelementToElement
(
current_element
);
// for every facet of the element
for
(
auto
&
current_facet
:
facets_to_element
)
{
if
(
current_facet
==
ElementNull
)
{
continue
;
}
if
(
current_facet
==
end_facet
)
{
facet_matched
=
true
;
}
auto
find_facet
=
std
::
find
(
facet_list
.
begin
(),
facet_list
.
end
(),
current_facet
);
// facet already listed or subfacet_connectivity is not in the
// connectivity of current_facet;
if
((
find_facet
!=
facet_list
.
end
())
or
not
hasElement
(
mesh_facets
.
getConnectivity
(
current_facet
),
subfacet_connectivity
))
{
continue
;
}
facet_list
.
push_back
(
current_facet
);
if
(
subfacet_list
!=
nullptr
)
{
const
Vector
<
Element
>
subfacets_of_facet
=
mesh_facets
.
getSubelementToElement
(
current_facet
);
/// check subfacets
for
(
const
auto
&
current_subfacet
:
subfacets_of_facet
)
{
if
(
current_subfacet
==
ElementNull
)
{
continue
;
}
auto
find_subfacet
=
std
::
find
(
subfacet_list
->
begin
(),
subfacet_list
->
end
(),
current_subfacet
);
if
((
find_subfacet
!=
subfacet_list
->
end
())
or
not
hasElement
(
mesh_facets
.
getConnectivity
(
current_subfacet
),
subfacet_connectivity
))
{
continue
;
}
subfacet_list
->
push_back
(
current_subfacet
);
}
}
/// consider opposing element
const
auto
&
elements_to_facet
=
mesh_facets
.
getElementToSubelement
(
current_facet
);
UInt
opposing
=
0
;
if
(
elements_to_facet
[
0
]
==
current_element
)
{
opposing
=
1
;
}
auto
&
opposing_element
=
elements_to_facet
[
opposing
];
/// skip null elements since they are on a boundary
if
(
opposing_element
==
ElementNull
)
{
continue
;
}
/// skip this element if already added
if
(
std
::
find
(
element_list
.
begin
(),
element_list
.
end
(),
opposing_element
)
!=
element_list
.
end
())
{
continue
;
}
/// only regular elements have to be checked
if
(
opposing_element
.
kind
()
==
_ek_regular
)
{
elements_to_check
.
push
(
opposing_element
);
}
element_list
.
push_back
(
opposing_element
);
AKANTU_DEBUG_ASSERT
(
hasElement
(
mesh
.
getConnectivity
(
opposing_element
),
subfacet_connectivity
),
"Subfacet doesn't belong to this element"
);
}
/// erased checked element from the list
elements_to_check
.
pop
();
}
return
facet_matched
;
}
/* -------------------------------------------------------------------------- */
void
CohesiveElementInserterHelper
::
updateCohesiveData
()
{
UInt
spatial_dimension
=
mesh
.
getSpatialDimension
();
bool
third_dimension
=
spatial_dimension
==
3
;
MeshAccessor
mesh_accessor
(
mesh
);
MeshAccessor
mesh_facets_accessor
(
mesh_facets
);
for
(
auto
ghost_type
:
ghost_types
)
{
for
(
auto
cohesive_type
:
mesh
.
elementTypes
(
_element_kind
=
_ek_cohesive
))
{
auto
nb_cohesive_elements
=
mesh
.
getNbElement
(
cohesive_type
,
ghost_type
);
nb_new_facets
(
cohesive_type
,
ghost_type
)
=
nb_cohesive_elements
;
mesh_facets_accessor
.
getSubelementToElement
(
cohesive_type
,
ghost_type
);
}
}
auto
&
facets_to_double
=
*
facets_to_double_by_dim
[
spatial_dimension
-
1
];
new_elements
.
reserve
(
new_elements
.
size
()
+
facets_to_double
.
size
());
std
::
array
<
Element
,
2
>
facets
;
auto
&
element_to_facet
=
mesh_facets_accessor
.
getElementToSubelement
();
auto
&
facets_to_cohesive_elements
=
mesh_facets_accessor
.
getSubelementToElement
();
for
(
auto
&&
facet_to_double
:
make_view
(
facets_to_double
,
2
))
{
auto
&
old_facet
=
facet_to_double
[
0
];
/// (in 3D cohesive elements connectivity is inverted)
facets
[
third_dimension
?
1
:
0
]
=
old_facet
;
facets
[
third_dimension
?
0
:
1
]
=
facet_to_double
[
1
];
auto
type_cohesive
=
FEEngine
::
getCohesiveElementType
(
old_facet
.
type
);
auto
&
facet_connectivities
=
mesh_facets
.
getConnectivity
(
old_facet
.
type
,
old_facet
.
ghost_type
);
auto
facet_connectivity_it
=
facet_connectivities
.
begin
(
facet_connectivities
.
getNbComponent
());
auto
cohesive_element
=
Element
{
type_cohesive
,
nb_new_facets
(
type_cohesive
,
old_facet
.
ghost_type
)
++
,
old_facet
.
ghost_type
};
auto
&
cohesives_connectivities
=
mesh_accessor
.
getConnectivity
(
type_cohesive
,
old_facet
.
ghost_type
);
Matrix
<
UInt
>
connectivity
(
facet_connectivities
.
getNbComponent
(),
2
);
Vector
<
Element
>
facets_to_cohesive_element
(
2
);
for
(
auto
s
:
arange
(
2
))
{
/// store doubled facets
facets_to_cohesive_element
[
s
]
=
facets
[
s
];
// update connectivities
connectivity
(
s
)
=
Vector
<
UInt
>
(
facet_connectivity_it
[
facets
[
s
].
element
]);
/// update element_to_facet vectors
element_to_facet
(
facets
[
s
],
0
)[
1
]
=
cohesive_element
;
}
cohesives_connectivities
.
push_back
(
connectivity
);
facets_to_cohesive_elements
(
type_cohesive
,
old_facet
.
ghost_type
)
.
push_back
(
facets_to_cohesive_element
);
/// add cohesive element to the element event list
new_elements
.
push_back
(
cohesive_element
);
}
}
/* -------------------------------------------------------------------------- */
void
CohesiveElementInserterHelper
::
doublePointFacet
()
{
UInt
spatial_dimension
=
mesh
.
getSpatialDimension
();
if
(
spatial_dimension
!=
1
)
{
return
;
}
auto
&
facets_to_double
=
*
facets_to_double_by_dim
[
spatial_dimension
-
1
];
auto
&
element_to_facet
=
mesh_facets
.
getElementToSubelement
();
auto
&
position
=
mesh
.
getNodes
();
MeshAccessor
mesh_accessor
(
mesh_facets
);
for
(
auto
ghost_type
:
ghost_types
)
{
for
(
auto
facet_type
:
nb_new_facets
.
elementTypes
(
spatial_dimension
-
1
,
ghost_type
,
_ek_regular
))
{
auto
nb_new_element
=
nb_new_facets
(
facet_type
,
ghost_type
);
auto
&
connectivities
=
mesh_accessor
.
getConnectivity
(
facet_type
,
ghost_type
);
connectivities
.
resize
(
connectivities
.
size
()
+
nb_new_element
);
}
}
for
(
auto
facet_to_double
:
make_view
(
facets_to_double
,
2
))
{
auto
&
old_facet
=
facet_to_double
[
0
];
auto
&
new_facet
=
facet_to_double
[
1
];
auto
element
=
element_to_facet
(
new_facet
)[
0
];
auto
&
facet_connectivities
=
mesh_accessor
.
getConnectivity
(
old_facet
.
type
,
old_facet
.
ghost_type
);
auto
old_node
=
facet_connectivities
(
old_facet
.
element
);
auto
new_node
=
position
.
size
();
/// update position
position
.
push_back
(
position
(
old_node
));
facet_connectivities
(
new_facet
.
element
)
=
new_node
;
Vector
<
UInt
>
segment_conectivity
=
mesh
.
getConnectivity
(
element
);
/// update facet connectivity
auto
it
=
std
::
find
(
segment_conectivity
.
begin
(),
segment_conectivity
.
end
(),
old_node
);
*
it
=
new_node
;
doubled_nodes
.
push_back
(
Vector
<
UInt
>
{
old_node
,
new_node
});
}
}
/* -------------------------------------------------------------------------- */
template
<
UInt
spatial_dimension
>
void
CohesiveElementInserterHelper
::
doubleSubfacet
()
{
if
(
spatial_dimension
==
1
)
{
return
;
}
std
::
vector
<
UInt
>
nodes_to_double
;
MeshAccessor
mesh_accessor
(
mesh_facets
);
auto
&
connectivities
=
mesh_accessor
.
getConnectivities
();
auto
&
facets_to_double
=
*
facets_to_double_by_dim
[
0
];
ElementTypeMap
<
Int
>
nb_new_elements
;
for
(
auto
&&
facet_to_double
:
make_view
(
facets_to_double
,
2
))
{
auto
&&
old_element
=
facet_to_double
[
0
];
nb_new_elements
(
old_element
.
type
,
old_element
.
ghost_type
)
=
0
;
}
for
(
auto
&&
facet_to_double
:
make_view
(
facets_to_double
,
2
))
{
auto
&&
old_element
=
facet_to_double
[
0
];
++
nb_new_elements
(
old_element
.
type
,
old_element
.
ghost_type
);
}
for
(
auto
ghost_type
:
ghost_types
)
{
for
(
auto
facet_type
:
nb_new_elements
.
elementTypes
(
0
,
ghost_type
,
_ek_regular
))
{
auto
&
connectivities
=
mesh_accessor
.
getConnectivity
(
facet_type
,
ghost_type
);
connectivities
.
resize
(
connectivities
.
size
()
+
nb_new_elements
(
facet_type
,
ghost_type
));
}
}
for
(
auto
&&
facet_to_double
:
make_view
(
facets_to_double
,
2
))
{
auto
&
old_facet
=
facet_to_double
(
0
);
// auto & new_facet = facet_to_double(1);
AKANTU_DEBUG_ASSERT
(
old_facet
.
type
==
_point_1
,
"Only _point_1 subfacet doubling is supported at the moment"
);
/// list nodes double
nodes_to_double
.
push_back
(
connectivities
(
old_facet
));
}
auto
old_nb_doubled_nodes
=
doubled_nodes
.
size
();
doubleNodes
(
nodes_to_double
);
auto
double_nodes_view
=
make_view
(
doubled_nodes
,
2
);
for
(
auto
&&
data
:
zip
(
make_view
(
facets_to_double
,
2
),
range
(
double_nodes_view
.
begin
()
+
old_nb_doubled_nodes
,
double_nodes_view
.
end
()),
arange
(
facets_to_double
.
size
())))
{
// auto & old_facet = std::get<0>(data)[0];
auto
&
new_facet
=
std
::
get
<
0
>
(
data
)[
1
];
auto
&
nodes
=
std
::
get
<
1
>
(
data
);
auto
old_node
=
nodes
(
0
);
auto
new_node
=
nodes
(
1
);
auto
f
=
std
::
get
<
2
>
(
data
);
/// add new nodes in connectivity
connectivities
(
new_facet
)
=
new_node
;
updateElementalConnectivity
(
mesh
,
old_node
,
new_node
,
elementsOfDimToElementsOfDim
(
2
,
0
)(
f
),
&
elementsOfDimToElementsOfDim
(
1
,
0
)(
f
));
updateElementalConnectivity
(
mesh_facets
,
old_node
,
new_node
,
elementsOfDimToElementsOfDim
(
1
,
0
)(
f
));
if
(
spatial_dimension
==
3
)
{
updateElementalConnectivity
(
mesh_facets
,
old_node
,
new_node
,
elementsOfDimToElementsOfDim
(
0
,
0
)(
f
));
}
}
updateSubelementToElement
(
0
,
spatial_dimension
==
2
);
updateElementToSubelement
(
0
,
spatial_dimension
==
2
);
}
}
// namespace akantu
Event Timeline
Log In to Comment