Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F78724778
mesh_igfem_spherical_growing_gel_tmpl.hh
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
Thu, Aug 22, 15:18
Size
19 KB
Mime Type
text/x-c++
Expires
Sat, Aug 24, 15:18 (2 d)
Engine
blob
Format
Raw Data
Handle
20084486
Attached To
rAKA akantu
mesh_igfem_spherical_growing_gel_tmpl.hh
View Options
/**
* @file mesh_igfem_spherical_growing_gel.cc
* @author Marco Vocialta <marco.vocialta@epfl.ch>
* @date Mon Sep 21 12:42:11 2015
*
* @brief Implementation of the functions of MeshIgfemSphericalGrowingGel
*
* @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 "mesh_utils.hh"
#include <algorithm>
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
template
<
UInt
dim
>
MeshIgfemSphericalGrowingGel
<
dim
>::
MeshIgfemSphericalGrowingGel
(
Mesh
&
mesh
)
:
mesh
(
mesh
),
nb_nodes_fem
(
mesh
.
getNbNodes
()),
nb_enriched_nodes
(
0
),
synchronizer
(
NULL
)
{
}
/* -------------------------------------------------------------------------- */
template
<
UInt
dim
>
void
MeshIgfemSphericalGrowingGel
<
dim
>::
init
()
{
nb_nodes_fem
=
mesh
.
getNbNodes
();
for
(
ghost_type_t
::
iterator
gt
=
ghost_type_t
::
begin
();
gt
!=
ghost_type_t
::
end
();
++
gt
)
{
GhostType
ghost_type
=
*
gt
;
Mesh
::
type_iterator
tit
=
mesh
.
firstType
(
dim
,
ghost_type
);
Mesh
::
type_iterator
tend
=
mesh
.
lastType
(
dim
,
ghost_type
);
for
(;
tit
!=
tend
;
++
tit
)
{
// loop to add corresponding IGFEM element types
if
(
*
tit
==
_triangle_3
){
mesh
.
addConnectivityType
(
_igfem_triangle_4
,
ghost_type
);
mesh
.
addConnectivityType
(
_igfem_triangle_5
,
ghost_type
);
}
else
AKANTU_DEBUG_ERROR
(
"Not ready for mesh type "
<<
*
tit
);
}
tit
=
mesh
.
firstType
(
dim
,
ghost_type
,
_ek_not_defined
);
tend
=
mesh
.
lastType
(
dim
,
ghost_type
,
_ek_not_defined
);
for
(;
tit
!=
tend
;
++
tit
)
{
AKANTU_BOOST_LIST_SWITCH
(
INSTANTIATOR
,
ELEMENT_LIST
,
*
tit
);
}
}
}
/* -------------------------------------------------------------------------- */
template
<
UInt
dim
>
void
MeshIgfemSphericalGrowingGel
<
dim
>::
computeMeshQueryListIntersectionPoint
(
const
std
::
list
<
SK
::
Sphere_3
>
&
query_list
)
{
/// store number of currently enriched nodes
this
->
nb_enriched_nodes
=
mesh
.
getNbNodes
()
-
nb_nodes_fem
;
UInt
nb_old_nodes
=
mesh
.
getNbNodes
();
for
(
ghost_type_t
::
iterator
gt
=
ghost_type_t
::
begin
();
gt
!=
ghost_type_t
::
end
();
++
gt
)
{
GhostType
ghost_type
=
*
gt
;
Mesh
::
type_iterator
iit
=
mesh
.
firstType
(
dim
,
ghost_type
,
_ek_not_defined
);
Mesh
::
type_iterator
iend
=
mesh
.
lastType
(
dim
,
ghost_type
,
_ek_not_defined
);
for
(;
iit
!=
iend
;
++
iit
)
{
MeshAbstractIntersector
<
SK
::
Sphere_3
>
&
intersector
=
*
intersectors
(
*
iit
,
ghost_type
);
intersector
.
constructData
(
ghost_type
);
intersector
.
computeMeshQueryListIntersectionPoint
(
query_list
,
nb_old_nodes
);
const
Array
<
Real
>
intersection_points_current_type
=
*
(
intersector
.
getIntersectionPoints
());
const
Array
<
UInt
>
&
new_node_per_elem
=
intersector
.
getNewNodePerElem
();
/// Send the new node event
UInt
new_points
=
intersection_points_current_type
.
getSize
();
StaticCommunicator
::
getStaticCommunicator
().
allReduce
(
&
new_points
,
1
,
_so_sum
);
if
(
new_points
>
0
)
{
Array
<
Real
>
&
nodes
=
this
->
mesh
.
getNodes
();
UInt
nb_node
=
nodes
.
getSize
();
Array
<
Real
>::
const_vector_iterator
intersec_p_it
=
intersection_points_current_type
.
begin
(
dim
);
NewIGFEMNodesEvent
new_nodes_event
;
for
(
UInt
in
=
0
;
in
<
intersection_points_current_type
.
getSize
();
++
in
,
++
intersec_p_it
)
{
nodes
.
push_back
(
*
intersec_p_it
);
new_nodes_event
.
getList
().
push_back
(
nb_node
+
in
);
}
new_nodes_event
.
setNewNodePerElem
(
new_node_per_elem
);
new_nodes_event
.
setType
(
*
iit
);
new_nodes_event
.
setGhostType
(
ghost_type
);
this
->
mesh
.
sendEvent
(
new_nodes_event
);
}
}
}
}
/* -------------------------------------------------------------------------- */
template
<
UInt
dim
>
void
MeshIgfemSphericalGrowingGel
<
dim
>::
removeAdditionalNodes
(){
AKANTU_DEBUG_IN
();
UInt
total_nodes
=
this
->
mesh
.
getNbNodes
();
UInt
nb_new_enriched_nodes
=
total_nodes
-
this
->
nb_enriched_nodes
-
this
->
nb_nodes_fem
;
UInt
old_total_nodes
=
this
->
nb_nodes_fem
+
this
->
nb_enriched_nodes
;
UInt
total_new_nodes
=
nb_new_enriched_nodes
;
StaticCommunicator
::
getStaticCommunicator
().
allReduce
(
&
total_new_nodes
,
1
,
_so_sum
);
if
(
total_new_nodes
==
0
)
return
;
RemovedNodesEvent
remove_nodes
(
this
->
mesh
);
Array
<
UInt
>
&
nodes_removed
=
remove_nodes
.
getList
();
Array
<
UInt
>
&
new_numbering
=
remove_nodes
.
getNewNumbering
();
for
(
UInt
nnod
=
0
;
nnod
<
this
->
nb_nodes_fem
;
++
nnod
){
new_numbering
(
nnod
)
=
nnod
;
}
for
(
UInt
nnod
=
nb_nodes_fem
;
nnod
<
old_total_nodes
;
++
nnod
){
new_numbering
(
nnod
)
=
UInt
(
-
1
)
;
nodes_removed
.
push_back
(
nnod
);
}
for
(
UInt
nnod
=
0
;
nnod
<
nb_new_enriched_nodes
;
++
nnod
){
new_numbering
(
nnod
+
old_total_nodes
)
=
this
->
nb_nodes_fem
+
nnod
;
}
if
(
nb_new_enriched_nodes
>
0
)
{
for
(
UInt
gt
=
_not_ghost
;
gt
<=
_ghost
;
++
gt
)
{
GhostType
ghost_type
=
(
GhostType
)
gt
;
Mesh
::
type_iterator
it
=
mesh
.
firstType
(
_all_dimensions
,
ghost_type
,
_ek_not_defined
);
Mesh
::
type_iterator
end
=
mesh
.
lastType
(
_all_dimensions
,
ghost_type
,
_ek_not_defined
);
for
(;
it
!=
end
;
++
it
)
{
ElementType
type
=
*
it
;
Array
<
UInt
>
&
connectivity_array
=
mesh
.
getConnectivity
(
type
,
ghost_type
);
UInt
nb_nodes_per_element
=
connectivity_array
.
getNbComponent
();
Array
<
UInt
>::
vector_iterator
conn_it
=
connectivity_array
.
begin
(
nb_nodes_per_element
);
Array
<
UInt
>::
vector_iterator
conn_end
=
connectivity_array
.
end
(
nb_nodes_per_element
);
for
(;
conn_it
!=
conn_end
;
++
conn_it
)
{
for
(
UInt
n
=
0
;
n
<
nb_nodes_per_element
;
++
n
)
{
(
*
conn_it
)(
n
)
=
new_numbering
((
*
conn_it
)(
n
));
}
}
}
}
}
this
->
mesh
.
sendEvent
(
remove_nodes
);
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
UInt
dim
>
void
MeshIgfemSphericalGrowingGel
<
dim
>::
buildIGFEMMesh
(){
AKANTU_DEBUG_IN
();
NewIGFEMElementsEvent
new_elements_event
;
UInt
total_new_elements
=
0
;
Array
<
Element
>
&
new_elements_list
=
new_elements_event
.
getList
();
Array
<
Element
>
&
old_elements_list
=
new_elements_event
.
getOldElementsList
();
RemovedElementsEvent
removed_elements_event
(
this
->
mesh
);
ChangedElementsEvent
changed_elements_event
(
this
->
mesh
);
for
(
ghost_type_t
::
iterator
gt
=
ghost_type_t
::
begin
();
gt
!=
ghost_type_t
::
end
();
++
gt
)
{
GhostType
ghost_type
=
*
gt
;
mesh
.
initElementTypeMapArray
(
removed_elements_event
.
getNewNumbering
(),
1
,
dim
,
ghost_type
,
false
,
_ek_not_defined
);
mesh
.
initElementTypeMapArray
(
changed_elements_event
.
getNewNumbering
(),
1
,
dim
,
ghost_type
,
false
,
_ek_not_defined
);
/// loop over all types in the mesh for a given ghost type
Mesh
::
type_iterator
iit
=
mesh
.
firstType
(
dim
,
ghost_type
,
_ek_not_defined
);
Mesh
::
type_iterator
iend
=
mesh
.
lastType
(
dim
,
ghost_type
,
_ek_not_defined
);
for
(;
iit
!=
iend
;
++
iit
)
{
ElementType
type
=
*
iit
;
MeshAbstractIntersector
<
SK
::
Sphere_3
>
&
intersector
=
*
intersectors
(
type
,
ghost_type
);
const
Array
<
UInt
>
&
new_node_per_elem
=
intersector
.
getNewNodePerElem
();
UInt
n_new_el
=
0
;
Array
<
UInt
>
&
connectivity
=
this
->
mesh
.
getConnectivity
(
type
,
ghost_type
);
/// get the connectivities of all types that that may transform
Array
<
UInt
>
&
connec_igfem_tri_4
=
this
->
mesh
.
getConnectivity
(
_igfem_triangle_4
,
ghost_type
);
Array
<
UInt
>
&
connec_igfem_tri_5
=
this
->
mesh
.
getConnectivity
(
_igfem_triangle_5
,
ghost_type
);
Array
<
UInt
>
&
connec_tri_3
=
this
->
mesh
.
getConnectivity
(
_triangle_3
,
ghost_type
);
/// create elements to store the newly generated elements
Element
el_tri_3
(
_triangle_3
,
0
,
ghost_type
,
_ek_regular
);
Element
el_igfem_tri_4
(
_igfem_triangle_4
,
0
,
ghost_type
,
_ek_igfem
);
Element
el_igfem_tri5
(
_igfem_triangle_5
,
0
,
ghost_type
,
_ek_igfem
);
Array
<
UInt
>
&
new_numbering
=
removed_elements_event
.
getNewNumbering
(
type
,
ghost_type
);
new_numbering
.
resize
(
connectivity
.
getSize
());
/// container for element to be removed
Element
old_element
(
type
,
0
,
ghost_type
,
Mesh
::
getKind
(
type
));
for
(
UInt
nel
=
0
;
nel
!=
new_node_per_elem
.
getSize
();
++
nel
)
{
/// a former IGFEM triangle is transformed into a regular triangle
if
(
(
type
!=
_triangle_3
)
&&
(
new_node_per_elem
(
nel
,
0
)
==
0
))
{
Vector
<
UInt
>
connec_new_elem
(
3
);
connec_new_elem
(
0
)
=
connectivity
(
nel
,
0
);
connec_new_elem
(
1
)
=
connectivity
(
nel
,
1
);
connec_new_elem
(
2
)
=
connectivity
(
nel
,
2
);
/// add the new element in the mesh
UInt
nb_triangles_3
=
connec_tri_3
.
getSize
();
connec_tri_3
.
push_back
(
connec_new_elem
);
el_tri_3
.
element
=
nb_triangles_3
;
new_elements_list
.
push_back
(
el_tri_3
);
/// the number of the old element in the mesh
old_element
.
element
=
nel
;
old_elements_list
.
push_back
(
old_element
);
new_numbering
(
nel
)
=
UInt
(
-
1
);
++
n_new_el
;
}
/// element is enriched
else
if
(
new_node_per_elem
(
nel
,
0
)
!=
0
)
{
switch
(
new_node_per_elem
(
nel
,
0
)){
/// new element is of type igfem_triangle_4
case
1
:
{
Vector
<
UInt
>
connec_new_elem
(
4
);
switch
(
new_node_per_elem
(
nel
,
2
)){
case
0
:
connec_new_elem
(
0
)
=
connectivity
(
nel
,
2
);
connec_new_elem
(
1
)
=
connectivity
(
nel
,
0
);
connec_new_elem
(
2
)
=
connectivity
(
nel
,
1
);
break
;
case
1
:
connec_new_elem
(
0
)
=
connectivity
(
nel
,
0
);
connec_new_elem
(
1
)
=
connectivity
(
nel
,
1
);
connec_new_elem
(
2
)
=
connectivity
(
nel
,
2
);
break
;
case
2
:
connec_new_elem
(
0
)
=
connectivity
(
nel
,
1
);
connec_new_elem
(
1
)
=
connectivity
(
nel
,
2
);
connec_new_elem
(
2
)
=
connectivity
(
nel
,
0
);
break
;
default
:
AKANTU_DEBUG_ERROR
(
"A triangle has only 3 segments not "
<<
new_node_per_elem
(
nel
,
2
));
break
;
}
connec_new_elem
(
3
)
=
new_node_per_elem
(
nel
,
1
);
UInt
nb_igfem_triangles_4
=
connec_igfem_tri_4
.
getSize
();
connec_igfem_tri_4
.
push_back
(
connec_new_elem
);
el_igfem_tri_4
.
element
=
nb_igfem_triangles_4
;
new_elements_list
.
push_back
(
el_igfem_tri_4
);
break
;
}
/// new element is of type igfem_triangle_5
case
2
:
{
Vector
<
UInt
>
connec_new_elem
(
5
);
if
(
(
new_node_per_elem
(
nel
,
2
)
==
0
)
&&
(
new_node_per_elem
(
nel
,
4
)
==
1
)
){
connec_new_elem
(
0
)
=
connectivity
(
nel
,
1
);
connec_new_elem
(
1
)
=
connectivity
(
nel
,
2
);
connec_new_elem
(
2
)
=
connectivity
(
nel
,
0
);
connec_new_elem
(
3
)
=
new_node_per_elem
(
nel
,
3
);
connec_new_elem
(
4
)
=
new_node_per_elem
(
nel
,
1
);
}
else
if
((
new_node_per_elem
(
nel
,
2
)
==
0
)
&&
(
new_node_per_elem
(
nel
,
4
)
==
2
)){
connec_new_elem
(
0
)
=
connectivity
(
nel
,
0
);
connec_new_elem
(
1
)
=
connectivity
(
nel
,
1
);
connec_new_elem
(
2
)
=
connectivity
(
nel
,
2
);
connec_new_elem
(
3
)
=
new_node_per_elem
(
nel
,
1
);
connec_new_elem
(
4
)
=
new_node_per_elem
(
nel
,
3
);
}
else
if
((
new_node_per_elem
(
nel
,
2
)
==
1
)
&&
(
new_node_per_elem
(
nel
,
4
)
==
2
)){
connec_new_elem
(
0
)
=
connectivity
(
nel
,
2
);
connec_new_elem
(
1
)
=
connectivity
(
nel
,
0
);
connec_new_elem
(
2
)
=
connectivity
(
nel
,
1
);
connec_new_elem
(
3
)
=
new_node_per_elem
(
nel
,
3
);
connec_new_elem
(
4
)
=
new_node_per_elem
(
nel
,
1
);
}
else
if
((
new_node_per_elem
(
nel
,
2
)
==
1
)
&&
(
new_node_per_elem
(
nel
,
4
)
==
0
)){
connec_new_elem
(
0
)
=
connectivity
(
nel
,
1
);
connec_new_elem
(
1
)
=
connectivity
(
nel
,
2
);
connec_new_elem
(
2
)
=
connectivity
(
nel
,
0
);
connec_new_elem
(
3
)
=
new_node_per_elem
(
nel
,
1
);
connec_new_elem
(
4
)
=
new_node_per_elem
(
nel
,
3
);
}
else
if
((
new_node_per_elem
(
nel
,
2
)
==
2
)
&&
(
new_node_per_elem
(
nel
,
4
)
==
0
)){
connec_new_elem
(
0
)
=
connectivity
(
nel
,
0
);
connec_new_elem
(
1
)
=
connectivity
(
nel
,
1
);
connec_new_elem
(
2
)
=
connectivity
(
nel
,
2
);
connec_new_elem
(
3
)
=
new_node_per_elem
(
nel
,
3
);
connec_new_elem
(
4
)
=
new_node_per_elem
(
nel
,
1
);
}
else
if
((
new_node_per_elem
(
nel
,
2
)
==
2
)
&&
(
new_node_per_elem
(
nel
,
4
)
==
1
)){
connec_new_elem
(
0
)
=
connectivity
(
nel
,
2
);
connec_new_elem
(
1
)
=
connectivity
(
nel
,
0
);
connec_new_elem
(
2
)
=
connectivity
(
nel
,
1
);
connec_new_elem
(
3
)
=
new_node_per_elem
(
nel
,
1
);
connec_new_elem
(
4
)
=
new_node_per_elem
(
nel
,
3
);
}
else
AKANTU_DEBUG_ERROR
(
"A triangle has only 3 segments (0 to 2) not "
<<
new_node_per_elem
(
nel
,
2
)
<<
"and"
<<
new_node_per_elem
(
nel
,
4
));
UInt
nb_igfem_triangles_5
=
connec_igfem_tri_5
.
getSize
();
connec_igfem_tri_5
.
push_back
(
connec_new_elem
);
el_igfem_tri5
.
element
=
nb_igfem_triangles_5
;
new_elements_list
.
push_back
(
el_igfem_tri5
);
break
;
}
default
:
AKANTU_DEBUG_ERROR
(
"IGFEM cannot add "
<<
new_node_per_elem
(
nel
,
0
)
<<
" nodes to triangles"
);
break
;
}
old_element
.
element
=
nel
;
old_elements_list
.
push_back
(
old_element
);
new_numbering
(
nel
)
=
UInt
(
-
1
);
++
n_new_el
;
}
}
total_new_elements
+=
n_new_el
;
}
}
/// update new numbering
for
(
ghost_type_t
::
iterator
gt
=
ghost_type_t
::
begin
();
gt
!=
ghost_type_t
::
end
();
++
gt
)
{
GhostType
ghost_type
=
*
gt
;
/// loop over all types in the mesh for a given ghost type
Mesh
::
type_iterator
iit
=
mesh
.
firstType
(
dim
,
ghost_type
,
_ek_not_defined
);
Mesh
::
type_iterator
iend
=
mesh
.
lastType
(
dim
,
ghost_type
,
_ek_not_defined
);
for
(;
iit
!=
iend
;
++
iit
)
{
ElementType
type
=
*
iit
;
Array
<
UInt
>
&
new_numbering
=
removed_elements_event
.
getNewNumbering
(
type
,
ghost_type
);
UInt
el_index
=
0
;
UInt
nb_element
=
this
->
mesh
.
getNbElement
(
type
,
ghost_type
);
new_numbering
.
resize
(
nb_element
);
for
(
UInt
e
=
0
;
e
<
nb_element
;
++
e
)
{
if
(
new_numbering
(
e
)
!=
UInt
(
-
1
))
{
new_numbering
(
e
)
=
el_index
;
++
el_index
;
}
}
changed_elements_event
.
getNewNumbering
(
type
,
ghost_type
).
copy
(
new_numbering
);
}
}
StaticCommunicator
::
getStaticCommunicator
().
allReduce
(
&
total_new_elements
,
1
,
_so_sum
);
if
(
total_new_elements
>
0
){
changed_elements_event
.
getListOld
().
copy
(
new_elements_event
.
getOldElementsList
());
changed_elements_event
.
getListNew
().
copy
(
new_elements_event
.
getList
());
this
->
mesh
.
sendEvent
(
changed_elements_event
);
this
->
mesh
.
sendEvent
(
new_elements_event
);
Array
<
Element
>
&
removed_list
=
removed_elements_event
.
getList
();
removed_list
.
copy
(
new_elements_event
.
getOldElementsList
());
this
->
mesh
.
sendEvent
(
removed_elements_event
);
}
removeAdditionalNodes
();
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
UInt
dim
>
void
MeshIgfemSphericalGrowingGel
<
dim
>::
buildSegmentConnectivityToNodeType
()
{
Mesh
mesh_facets
(
mesh
.
initMeshFacets
());
MeshUtils
::
buildSegmentToNodeType
(
mesh
,
mesh_facets
,
synchronizer
);
// only the ghost elements are considered
for
(
UInt
g
=
_not_ghost
;
g
<=
_ghost
;
++
g
)
{
GhostType
ghost_type
=
(
GhostType
)
g
;
Mesh
::
type_iterator
it
=
mesh_facets
.
firstType
(
1
,
ghost_type
);
Mesh
::
type_iterator
end
=
mesh_facets
.
lastType
(
1
,
ghost_type
);
for
(;
it
!=
end
;
++
it
)
{
ElementType
type
=
*
it
;
const
Array
<
Int
>
&
segment_to_nodetype
=
mesh_facets
.
getData
<
Int
>
(
"segment_to_nodetype"
,
type
,
ghost_type
);
const
Array
<
UInt
>
&
segment_connectivity
=
mesh_facets
.
getConnectivity
(
type
,
ghost_type
);
// looping over all the segments
Array
<
UInt
>::
const_vector_iterator
conn_it
=
segment_connectivity
.
begin
(
segment_connectivity
.
getNbComponent
());
Array
<
UInt
>::
const_vector_iterator
conn_end
=
segment_connectivity
.
end
(
segment_connectivity
.
getNbComponent
());
UInt
seg_index
=
0
;
UInt
n
[
2
];
for
(;
conn_it
!=
conn_end
;
++
conn_it
,
++
seg_index
)
{
Int
seg_type
=
segment_to_nodetype
(
seg_index
);
n
[
0
]
=
(
*
conn_it
)(
0
);
n
[
1
]
=
(
*
conn_it
)(
1
);
if
(
(
mesh
.
isMasterNode
(
n
[
0
])
||
mesh
.
isSlaveNode
(
n
[
0
]))
&&
(
mesh
.
isMasterNode
(
n
[
1
])
||
mesh
.
isSlaveNode
(
n
[
1
])))
{
std
::
sort
(
n
,
n
+
2
);
segment_conn_to_node_type
[
std
::
make_pair
(
n
[
0
],
n
[
1
])]
=
seg_type
;
}
}
}
}
}
/* -------------------------------------------------------------------------- */
template
<
UInt
dim
>
void
MeshIgfemSphericalGrowingGel
<
dim
>::
updateNodeType
(
const
Array
<
UInt
>
&
nodes_list
,
const
Array
<
UInt
>
&
new_node_per_elem
,
ElementType
type
,
GhostType
ghost_type
)
{
Array
<
Int
>
&
nodes_type
=
mesh
.
getNodesType
();
UInt
old_nodes
=
nodes_type
.
getSize
();
UInt
new_nodes
=
nodes_list
.
getSize
();
// exit this function if the simulation in run in serial
if
(
old_nodes
==
0
||
new_nodes
==
0
)
return
;
nodes_type
.
resize
(
old_nodes
+
new_nodes
);
Array
<
bool
>
checked_node
(
new_nodes
,
1
,
false
);
UInt
nb_prim_by_el
=
0
;
if
(
(
type
==
_triangle_3
)
||
(
type
==
_igfem_triangle_4
)
||
(
type
==
_igfem_triangle_5
)
){
nb_prim_by_el
=
3
;
}
else
{
AKANTU_DEBUG_ERROR
(
"Not ready for mesh type "
<<
type
);
}
// determine the node type for the local, master and slave nodes
const
Array
<
UInt
>
&
connectivity
=
mesh
.
getConnectivity
(
type
,
ghost_type
);
unordered_map
<
std
::
pair
<
UInt
,
UInt
>
,
Int
>::
type
::
iterator
seg_type_it
;
unordered_map
<
std
::
pair
<
UInt
,
UInt
>
,
Int
>::
type
::
iterator
seg_type_end
=
segment_conn_to_node_type
.
end
();
for
(
UInt
el
=
0
;
el
<
new_node_per_elem
.
getSize
();
++
el
)
{
UInt
nb_nodes
=
new_node_per_elem
(
el
,
0
);
for
(
UInt
n
=
0
;
n
<
nb_nodes
;
++
n
)
{
UInt
node_index
=
new_node_per_elem
(
el
,
2
*
n
+
1
);
if
(
node_index
<
old_nodes
||
checked_node
(
node_index
-
old_nodes
))
continue
;
// get the elemental connectivity of the segment associated to the node
UInt
segment_index
=
new_node_per_elem
(
el
,
2
*
n
+
2
);
UInt
extreme_nodes
[
2
];
extreme_nodes
[
0
]
=
segment_index
;
extreme_nodes
[
1
]
=
segment_index
+
1
;
if
(
extreme_nodes
[
1
]
==
nb_prim_by_el
)
extreme_nodes
[
1
]
=
0
;
// get the connectivity of the segment
extreme_nodes
[
0
]
=
connectivity
(
el
,
extreme_nodes
[
0
]);
extreme_nodes
[
1
]
=
connectivity
(
el
,
extreme_nodes
[
1
]);
// if one extreme nodes is pure ghost, then also the new node is pure ghost
if
(
mesh
.
isPureGhostNode
(
extreme_nodes
[
0
])
||
mesh
.
isPureGhostNode
(
extreme_nodes
[
1
]))
nodes_type
(
node_index
)
=
-
3
;
// if on of the two extreme nodes is local, then also the new node is local
else
if
(
mesh
.
isLocalNode
(
extreme_nodes
[
0
])
||
mesh
.
isLocalNode
(
extreme_nodes
[
1
]))
nodes_type
(
node_index
)
=
-
1
;
// otherwise use the value stored in the map
else
{
std
::
sort
(
extreme_nodes
,
extreme_nodes
+
2
);
seg_type_it
=
segment_conn_to_node_type
.
find
(
std
::
make_pair
(
extreme_nodes
[
0
],
extreme_nodes
[
1
]));
AKANTU_DEBUG_ASSERT
(
seg_type_it
!=
seg_type_end
,
"Segment not found"
);
nodes_type
(
node_index
)
=
seg_type_it
->
second
;
}
checked_node
(
node_index
-
old_nodes
)
=
true
;
}
}
AKANTU_DEBUG_ASSERT
(
std
::
accumulate
(
checked_node
.
begin
(),
checked_node
.
end
(),
0
)
==
Int
(
checked_node
.
getSize
()),
"Not all new nodes were assigned a node type"
);
}
__END_AKANTU__
Event Timeline
Log In to Comment