Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F63405783
material_non_local_inline_impl.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
Sun, May 19, 21:57
Size
51 KB
Mime Type
text/x-c
Expires
Tue, May 21, 21:57 (2 d)
Engine
blob
Format
Raw Data
Handle
17769083
Attached To
rAKA akantu
material_non_local_inline_impl.cc
View Options
/**
* @file material_non_local_inline_impl.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date Wed Aug 31 11:09:48 2011
*
* @brief Non-local inline implementation
*
* @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/>.
*
*/
/* -------------------------------------------------------------------------- */
__END_AKANTU__
/* -------------------------------------------------------------------------- */
#include "aka_types.hh"
#include "grid_synchronizer.hh"
#include "synchronizer_registry.hh"
#include "integrator.hh"
/* -------------------------------------------------------------------------- */
#include <fstream>
/* -------------------------------------------------------------------------- */
#if defined(AKANTU_DEBUG_TOOLS)
# include "aka_debug_tools.hh"
#endif
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
template
<
UInt
DIM
,
template
<
UInt
>
class
WeightFunction
>
MaterialNonLocal
<
DIM
,
WeightFunction
>::
MaterialNonLocal
(
SolidMechanicsModel
&
model
,
const
ID
&
id
)
:
Material
(
model
,
id
),
weight_func
(
NULL
),
spatial_grid
(
NULL
),
compute_stress_calls
(
0
),
is_creating_grid
(
false
),
grid_synchronizer
(
NULL
)
{
AKANTU_DEBUG_IN
();
this
->
is_non_local
=
true
;
this
->
weight_func
=
new
WeightFunction
<
DIM
>
(
*
this
);
this
->
registerSubSection
(
_st_non_local
,
"weight_function"
,
*
weight_func
);
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
UInt
spatial_dimension
,
template
<
UInt
>
class
WeightFunction
>
MaterialNonLocal
<
spatial_dimension
,
WeightFunction
>::~
MaterialNonLocal
()
{
AKANTU_DEBUG_IN
();
delete
spatial_grid
;
delete
weight_func
;
delete
grid_synchronizer
;
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
UInt
spatial_dimension
,
template
<
UInt
>
class
WeightFunction
>
void
MaterialNonLocal
<
spatial_dimension
,
WeightFunction
>::
initMaterial
()
{
AKANTU_DEBUG_IN
();
// Material::initMaterial();
Mesh
&
mesh
=
this
->
model
->
getFEEngine
().
getMesh
();
InternalField
<
Real
>
quadrature_points_coordinates
(
"quadrature_points_coordinates_tmp_nl"
,
*
this
);
quadrature_points_coordinates
.
initialize
(
spatial_dimension
);
ElementTypeMap
<
UInt
>
nb_ghost_protected
;
Mesh
::
type_iterator
it
=
mesh
.
firstType
(
spatial_dimension
,
_ghost
);
Mesh
::
type_iterator
last_type
=
mesh
.
lastType
(
spatial_dimension
,
_ghost
);
for
(;
it
!=
last_type
;
++
it
)
nb_ghost_protected
(
mesh
.
getNbElement
(
*
it
,
_ghost
),
*
it
,
_ghost
);
AKANTU_DEBUG_INFO
(
"Creating cell list"
);
createCellList
(
quadrature_points_coordinates
);
AKANTU_DEBUG_INFO
(
"Creating pairs"
);
updatePairList
(
quadrature_points_coordinates
);
#if not defined(AKANTU_NDEBUG)
if
(
AKANTU_DEBUG_TEST
(
dblDump
))
neighbourhoodStatistics
(
"material_non_local.stats"
);
#endif
AKANTU_DEBUG_INFO
(
"Cleaning extra ghosts"
);
cleanupExtraGhostElement
(
nb_ghost_protected
);
AKANTU_DEBUG_INFO
(
"Computing weights"
);
weight_func
->
setRadius
(
weight_func
->
getRadius
());
weight_func
->
init
();
computeWeights
(
quadrature_points_coordinates
);
// >>>>>> DEBUG CODE >>>>>> //
#if defined(AKANTU_DEBUG_TOOLS) && defined(AKANTU_CORE_CXX11)
debug
::
element_manager
.
print
(
debug
::
_dm_material_non_local
,
[
this
,
&
quadrature_points_coordinates
](
const
Element
&
el
)
->
std
::
string
{
std
::
stringstream
out
;
FEEngine
&
fem
=
this
->
model
->
getFEEngine
();
Mesh
mesh
(
spatial_dimension
,
getID
()
+
"mesh_tmp"
);
mesh
.
addConnectivityType
(
_segment_2
);
Array
<
UInt
>
&
connectivity
=
const_cast
<
Array
<
UInt
>
&>
(
mesh
.
getConnectivity
(
_segment_2
));
Array
<
Real
>
&
nodes
=
const_cast
<
Array
<
Real
>
&>
(
mesh
.
getNodes
());
Array
<
bool
>
g_notg
(
0
,
1
);
Array
<
Real
>
dist
(
0
,
1
);
Array
<
Real
>
weight
(
0
,
1
);
Array
<
Real
>
damage
(
0
,
1
);
Array
<
Real
>
jac
(
0
,
1
);
UInt
nb_quad_per_elem
=
this
->
model
->
getFEEngine
().
getNbQuadraturePoints
(
el
.
type
,
el
.
ghost_type
);
Array
<
Real
>::
const_vector_iterator
quad_it
=
quadrature_points_coordinates
(
el
.
type
,
el
.
ghost_type
).
begin
(
spatial_dimension
);
std
::
map
<
Vector
<
Real
>
,
UInt
>
numbering
;
UInt
counter
=
0
;
const
Vector
<
Real
>
quad
=
quad_it
[
el
.
element
*
nb_quad_per_elem
];
numbering
[
el
.
element
*
nb_quad_per_elem
]
=
counter
++
;
nodes
.
push_back
(
quad
);
g_notg
.
push_back
(
el
.
ghost_type
==
_ghost
);
dist
.
push_back
(
quad
.
distance
(
quad
));
weight
.
push_back
(
0.
);
damage
.
push_back
((
getArray
(
"damage"
,
el
.
type
,
el
.
ghost_type
)(
el
.
element
*
fem
.
getNbQuadraturePoints
(
el
.
type
,
el
.
ghost_type
))));
jac
.
push_back
(
fem
.
getIntegratorInterface
().
getJacobians
(
el
.
type
,
el
.
ghost_type
)(
el
.
element
*
fem
.
getNbQuadraturePoints
(
el
.
type
,
el
.
ghost_type
)));
Vector
<
UInt
>
conn
(
2
);
conn
(
0
)
=
0
;
bool
found
=
false
;
GhostType
ghost_type1
=
_not_ghost
;
for
(
ghost_type_t
::
iterator
git
=
ghost_type_t
::
begin
();
git
!=
ghost_type_t
::
end
();
++
git
)
{
GhostType
ghost_type2
=
*
git
;
UInt
existing_pairs_num
=
ghost_type2
-
_not_ghost
;
pair_type
::
iterator
first_pair_types
=
existing_pairs
[
existing_pairs_num
].
begin
();
pair_type
::
iterator
last_pair_types
=
existing_pairs
[
existing_pairs_num
].
end
();
// Compute the weights
for
(;
first_pair_types
!=
last_pair_types
;
++
first_pair_types
)
{
ElementType
type1
=
first_pair_types
->
first
;
ElementType
type2
=
first_pair_types
->
second
;
const
Array
<
UInt
>
&
elem_mat_1
=
element_filter
(
type1
,
ghost_type1
);
const
Array
<
UInt
>
&
elem_mat_2
=
element_filter
(
type2
,
ghost_type2
);
const
Array
<
UInt
>
&
pairs
=
pair_list
(
type1
,
ghost_type1
)(
type2
,
ghost_type2
);
const
Array
<
Real
>
&
weights
=
pair_weight
(
type1
,
ghost_type1
)(
type2
,
ghost_type2
);
UInt
nb_quad1
=
fem
.
getNbQuadraturePoints
(
type1
,
ghost_type1
);
UInt
nb_quad2
=
fem
.
getNbQuadraturePoints
(
type2
,
ghost_type2
);
Array
<
UInt
>::
const_iterator
<
Vector
<
UInt
>
>
first_pair
=
pairs
.
begin
(
2
);
Array
<
UInt
>::
const_iterator
<
Vector
<
UInt
>
>
last_pair
=
pairs
.
end
(
2
);
Array
<
Real
>::
const_vector_iterator
pair_w
=
weights
.
begin
(
2
);
for
(;
first_pair
!=
last_pair
;
++
first_pair
,
++
pair_w
)
{
UInt
_q1
=
(
*
first_pair
)(
0
);
UInt
_q2
=
(
*
first_pair
)(
1
);
QuadraturePoint
q1
(
type1
,
elem_mat_1
(
_q1
/
nb_quad1
),
_q1
%
nb_quad1
,
ghost_type1
);
QuadraturePoint
q2
(
type2
,
elem_mat_2
(
_q2
/
nb_quad2
),
_q2
%
nb_quad2
,
ghost_type2
);
if
(
el
==
(
Element
)
q1
||
el
==
(
Element
)
q2
)
{
QuadraturePoint
q
;
Real
w1
;
UInt
ele
;
// Real w2;
if
(
el
==
(
Element
)
q1
)
{
q
=
q2
;
w1
=
(
*
pair_w
)(
0
);
ele
=
q
.
element
*
nb_quad2
;
// w2 = (*pair_w)(1);
}
else
{
q
=
q1
;
w1
=
(
*
pair_w
)(
1
);
ele
=
q
.
element
*
nb_quad1
;
// w2 = (*pair_w)(0);
}
found
=
true
;
fem
.
getNbQuadraturePoints
(
q
.
type
,
q
.
ghost_type
);
quad_it
=
quadrature_points_coordinates
(
q
.
type
,
q
.
ghost_type
).
begin
(
spatial_dimension
);
const
Vector
<
Real
>
&
quad_coord
=
quad_it
[
q
.
element
*
nb_quad_per_elem
];
std
::
map
<
Vector
<
Real
>
,
UInt
>::
iterator
nit
=
numbering
.
find
(
quad_coord
);
if
(
nit
!=
numbering
.
end
())
{
conn
(
1
)
=
nit
->
second
;
if
(
conn
(
1
)
==
0
)
{
weight
(
0
)
=
w1
;
}
}
else
{
conn
(
1
)
=
counter
++
;
numbering
[
quad_coord
]
=
conn
(
1
);
nodes
.
push_back
(
quad_coord
);
g_notg
.
push_back
(
q
.
ghost_type
==
_ghost
);
dist
.
push_back
(
quad
.
distance
(
quad_coord
));
weight
.
push_back
(
w1
);
damage
.
push_back
((
getArray
(
"damage"
,
q
.
type
,
q
.
ghost_type
)(
q
.
global_num
)));
jac
.
push_back
(
fem
.
getIntegratorInterface
().
getJacobians
(
q
.
type
,
q
.
ghost_type
)(
ele
));
}
connectivity
.
push_back
(
conn
);
}
}
}
}
if
(
found
)
{
std
::
stringstream
sstr
;
sstr
<<
"neigh_mesh"
<<
el
;
out
<<
"./neighbors"
<<
sstr
.
str
()
<<
" dumped!"
;
DumperParaview
dumper
(
sstr
.
str
(),
"./neighbors"
,
false
);
dumper
.
registerMesh
(
mesh
);
dumper
.
registerField
(
"ghost"
,
new
DumperIOHelper
::
NodalField
<
bool
>
(
g_notg
));
dumper
.
registerField
(
"distance"
,
new
DumperIOHelper
::
NodalField
<
Real
>
(
dist
));
dumper
.
registerField
(
"weight"
,
new
DumperIOHelper
::
NodalField
<
Real
>
(
weight
));
dumper
.
registerField
(
"damage"
,
new
DumperIOHelper
::
NodalField
<
Real
>
(
damage
));
dumper
.
registerField
(
"jacobian"
,
new
DumperIOHelper
::
NodalField
<
Real
>
(
jac
));
dumper
.
dump
();
}
return
out
.
str
();
});
#endif
// <<<<<< DEBUG CODE <<<<<< //
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
UInt
spatial_dimension
,
template
<
UInt
>
class
WeightFunction
>
void
MaterialNonLocal
<
spatial_dimension
,
WeightFunction
>::
cleanupExtraGhostElement
(
const
ElementTypeMap
<
UInt
>
&
nb_ghost_protected
)
{
AKANTU_DEBUG_IN
();
// Create list of element to keep
std
::
set
<
Element
>
relevant_ghost_element
;
pair_type
::
const_iterator
first_pair_types
=
existing_pairs
[
1
].
begin
();
pair_type
::
const_iterator
last_pair_types
=
existing_pairs
[
1
].
end
();
for
(;
first_pair_types
!=
last_pair_types
;
++
first_pair_types
)
{
ElementType
type2
=
first_pair_types
->
second
;
GhostType
ghost_type2
=
_ghost
;
UInt
nb_quad2
=
this
->
model
->
getFEEngine
().
getNbQuadraturePoints
(
type2
);
Array
<
UInt
>
&
elem_filter
=
element_filter
(
type2
,
ghost_type2
);
const
Array
<
UInt
>
&
pairs
=
pair_list
(
first_pair_types
->
first
,
_not_ghost
)(
first_pair_types
->
second
,
ghost_type2
);
Array
<
UInt
>::
const_iterator
<
Vector
<
UInt
>
>
first_pair
=
pairs
.
begin
(
2
);
Array
<
UInt
>::
const_iterator
<
Vector
<
UInt
>
>
last_pair
=
pairs
.
end
(
2
);
for
(;
first_pair
!=
last_pair
;
++
first_pair
)
{
UInt
_q2
=
(
*
first_pair
)(
1
);
QuadraturePoint
q2
(
type2
,
elem_filter
(
_q2
/
nb_quad2
),
_q2
%
nb_quad2
,
ghost_type2
);
relevant_ghost_element
.
insert
(
q2
);
}
}
// Create list of element to remove and new numbering for element to keep
Mesh
&
mesh
=
this
->
model
->
getFEEngine
().
getMesh
();
std
::
set
<
Element
>
ghost_to_erase
;
Mesh
::
type_iterator
it
=
mesh
.
firstType
(
spatial_dimension
,
_ghost
);
Mesh
::
type_iterator
last_type
=
mesh
.
lastType
(
spatial_dimension
,
_ghost
);
RemovedElementsEvent
remove_elem
(
mesh
);
Element
element
;
element
.
ghost_type
=
_ghost
;
for
(;
it
!=
last_type
;
++
it
)
{
element
.
type
=
*
it
;
UInt
nb_ghost_elem
=
mesh
.
getNbElement
(
*
it
,
_ghost
);
UInt
nb_ghost_elem_protected
=
0
;
try
{
nb_ghost_elem_protected
=
nb_ghost_protected
(
*
it
,
_ghost
);
}
catch
(...)
{}
if
(
!
remove_elem
.
getNewNumbering
().
exists
(
*
it
,
_ghost
))
remove_elem
.
getNewNumbering
().
alloc
(
nb_ghost_elem
,
1
,
*
it
,
_ghost
);
else
remove_elem
.
getNewNumbering
(
*
it
,
_ghost
).
resize
(
nb_ghost_elem
);
Array
<
UInt
>
&
elem_filter
=
element_filter
(
*
it
,
_ghost
);
Array
<
UInt
>
&
new_numbering
=
remove_elem
.
getNewNumbering
(
*
it
,
_ghost
);
UInt
ng
=
0
;
for
(
UInt
g
=
0
;
g
<
nb_ghost_elem
;
++
g
)
{
element
.
element
=
elem_filter
(
g
);
if
(
element
.
element
>=
nb_ghost_elem_protected
&&
(
std
::
find
(
relevant_ghost_element
.
begin
(),
relevant_ghost_element
.
end
(),
element
)
==
relevant_ghost_element
.
end
()))
{
ghost_to_erase
.
insert
(
element
);
remove_elem
.
getList
().
push_back
(
element
);
new_numbering
(
g
)
=
UInt
(
-
1
);
}
else
{
new_numbering
(
g
)
=
ng
;
++
ng
;
}
}
}
mesh
.
sendEvent
(
remove_elem
);
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
UInt
spatial_dimension
,
template
<
UInt
>
class
WeightFunction
>
void
MaterialNonLocal
<
spatial_dimension
,
WeightFunction
>::
createCellList
(
ElementTypeMapArray
<
Real
>
&
quadrature_points_coordinates
)
{
AKANTU_DEBUG_IN
();
const
Real
safety_factor
=
1.2
;
// for the cell grid spacing
Mesh
&
mesh
=
this
->
model
->
getFEEngine
().
getMesh
();
mesh
.
computeBoundingBox
();
Real
lower_bounds
[
spatial_dimension
];
Real
upper_bounds
[
spatial_dimension
];
mesh
.
getLocalLowerBounds
(
lower_bounds
);
mesh
.
getLocalUpperBounds
(
upper_bounds
);
Vector
<
Real
>
spacing
(
spatial_dimension
,
weight_func
->
getRadius
()
*
safety_factor
);
Vector
<
Real
>
center
(
spatial_dimension
);
for
(
UInt
i
=
0
;
i
<
spatial_dimension
;
++
i
)
{
center
(
i
)
=
(
upper_bounds
[
i
]
+
lower_bounds
[
i
])
/
2.
;
}
spatial_grid
=
new
SpatialGrid
<
QuadraturePoint
>
(
spatial_dimension
,
spacing
,
center
);
computeQuadraturePointsCoordinates
(
quadrature_points_coordinates
,
_not_ghost
);
fillCellList
(
quadrature_points_coordinates
,
_not_ghost
);
is_creating_grid
=
true
;
SynchronizerRegistry
&
synch_registry
=
this
->
model
->
getSynchronizerRegistry
();
std
::
stringstream
sstr
;
sstr
<<
getID
()
<<
":grid_synchronizer"
;
grid_synchronizer
=
GridSynchronizer
::
createGridSynchronizer
(
mesh
,
*
spatial_grid
,
sstr
.
str
());
synch_registry
.
registerSynchronizer
(
*
grid_synchronizer
,
_gst_mnl_for_average
);
synch_registry
.
registerSynchronizer
(
*
grid_synchronizer
,
_gst_mnl_weight
);
is_creating_grid
=
false
;
#if not defined(AKANTU_NDEBUG)
Mesh
*
mesh_tmp
=
NULL
;
if
(
AKANTU_DEBUG_TEST
(
dblDump
))
{
mesh_tmp
=
new
Mesh
(
spatial_dimension
,
"mnl_grid"
);
spatial_grid
->
saveAsMesh
(
*
mesh_tmp
);
std
::
stringstream
sstr_grid
;
StaticCommunicator
&
comm
=
StaticCommunicator
::
getStaticCommunicator
();
Int
prank
=
comm
.
whoAmI
();
sstr_grid
<<
"material_non_local_grid_"
<<
std
::
setfill
(
'0'
)
<<
std
::
setw
(
4
)
<<
prank
<<
".msh"
;
mesh_tmp
->
write
(
sstr_grid
.
str
());
delete
mesh_tmp
;
}
#endif
this
->
computeQuadraturePointsCoordinates
(
quadrature_points_coordinates
,
_ghost
);
fillCellList
(
quadrature_points_coordinates
,
_ghost
);
#if not defined(AKANTU_NDEBUG)
if
(
AKANTU_DEBUG_TEST
(
dblDump
))
{
mesh_tmp
=
new
Mesh
(
spatial_dimension
,
"mnl_grid"
);
spatial_grid
->
saveAsMesh
(
*
mesh_tmp
);
std
::
stringstream
sstr_grid
;
StaticCommunicator
&
comm
=
StaticCommunicator
::
getStaticCommunicator
();
Int
prank
=
comm
.
whoAmI
();
sstr_grid
.
str
(
std
::
string
());
sstr_grid
<<
"material_non_local_grid_ghost_"
<<
std
::
setfill
(
'0'
)
<<
std
::
setw
(
4
)
<<
prank
<<
".msh"
;
mesh_tmp
->
write
(
sstr_grid
.
str
());
delete
mesh_tmp
;
}
#endif
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
UInt
spatial_dimension
,
template
<
UInt
>
class
WeightFunction
>
void
MaterialNonLocal
<
spatial_dimension
,
WeightFunction
>::
fillCellList
(
const
ElementTypeMapArray
<
Real
>
&
quadrature_points_coordinates
,
const
GhostType
&
ghost_type
)
{
Mesh
&
mesh
=
this
->
model
->
getFEEngine
().
getMesh
();
QuadraturePoint
q
;
q
.
ghost_type
=
ghost_type
;
Mesh
::
type_iterator
it
=
mesh
.
firstType
(
spatial_dimension
,
ghost_type
);
Mesh
::
type_iterator
last_type
=
mesh
.
lastType
(
spatial_dimension
,
ghost_type
);
for
(;
it
!=
last_type
;
++
it
)
{
Array
<
UInt
>
&
elem_filter
=
element_filter
(
*
it
,
ghost_type
);
UInt
nb_element
=
elem_filter
.
getSize
();
UInt
nb_quad
=
this
->
model
->
getFEEngine
().
getNbQuadraturePoints
(
*
it
,
ghost_type
);
const
Array
<
Real
>
&
quads
=
quadrature_points_coordinates
(
*
it
,
ghost_type
);
q
.
type
=
*
it
;
Array
<
Real
>::
const_vector_iterator
quad
=
quads
.
begin
(
spatial_dimension
);
UInt
*
elem
=
elem_filter
.
storage
();
for
(
UInt
e
=
0
;
e
<
nb_element
;
++
e
)
{
q
.
element
=
*
elem
;
for
(
UInt
nq
=
0
;
nq
<
nb_quad
;
++
nq
)
{
q
.
num_point
=
nq
;
//q.setPosition(*quad);
spatial_grid
->
insert
(
q
,
*
quad
);
++
quad
;
}
++
elem
;
}
}
}
/* -------------------------------------------------------------------------- */
template
<
UInt
spatial_dimension
,
template
<
UInt
>
class
WeightFunction
>
void
MaterialNonLocal
<
spatial_dimension
,
WeightFunction
>::
updatePairList
(
const
ElementTypeMapArray
<
Real
>
&
quadrature_points_coordinates
)
{
AKANTU_DEBUG_IN
();
Mesh
&
mesh
=
this
->
model
->
getFEEngine
().
getMesh
();
GhostType
ghost_type
=
_not_ghost
;
QuadraturePoint
quad_point
;
quad_point
.
ghost_type
=
ghost_type
;
// generate the pair of neighbor depending of the cell_list
Mesh
::
type_iterator
it
=
mesh
.
firstType
(
spatial_dimension
,
ghost_type
);
Mesh
::
type_iterator
last_type
=
mesh
.
lastType
(
spatial_dimension
,
ghost_type
);
for
(;
it
!=
last_type
;
++
it
)
{
// Preparing datas
const
Array
<
Real
>
&
quads
=
quadrature_points_coordinates
(
*
it
,
ghost_type
);
Array
<
Real
>::
const_vector_iterator
first_quad
=
quads
.
begin
(
spatial_dimension
);
Array
<
Real
>::
const_vector_iterator
last_quad
=
quads
.
end
(
spatial_dimension
);
ElementTypeMapArray
<
UInt
>
&
pairs
=
pair_list
(
ElementTypeMapArray
<
UInt
>
(
"pairs"
,
getID
(),
memory_id
),
*
it
,
ghost_type
);
ElementType
current_element_type
=
_not_defined
;
GhostType
current_ghost_type
=
_casper
;
UInt
existing_pairs_num
=
0
;
Array
<
UInt
>
*
neighbors
=
NULL
;
Array
<
UInt
>::
const_iterator
<
Vector
<
UInt
>
>
element_index_material_it
;
Array
<
Real
>::
const_vector_iterator
quad_coord_it
;
UInt
my_num_quad
=
0
;
quad_point
.
type
=
*
it
;
// loop over quad points
for
(;
first_quad
!=
last_quad
;
++
first_quad
,
++
my_num_quad
)
{
SpatialGrid
<
QuadraturePoint
>::
CellID
cell_id
=
spatial_grid
->
getCellID
(
*
first_quad
);
SpatialGrid
<
QuadraturePoint
>::
neighbor_cells_iterator
first_neigh_cell
=
spatial_grid
->
beginNeighborCells
(
cell_id
);
SpatialGrid
<
QuadraturePoint
>::
neighbor_cells_iterator
last_neigh_cell
=
spatial_grid
->
endNeighborCells
(
cell_id
);
quad_point
.
element
=
my_num_quad
/
this
->
model
->
getFEEngine
().
getNbQuadraturePoints
(
*
it
,
quad_point
.
ghost_type
);
// loop over neighbors cells of the one containing the current quadrature
// point
for
(;
first_neigh_cell
!=
last_neigh_cell
;
++
first_neigh_cell
)
{
SpatialGrid
<
QuadraturePoint
>::
Cell
::
iterator
first_neigh_quad
=
spatial_grid
->
beginCell
(
*
first_neigh_cell
);
SpatialGrid
<
QuadraturePoint
>::
Cell
::
iterator
last_neigh_quad
=
spatial_grid
->
endCell
(
*
first_neigh_cell
);
// loop over the quadrature point in the current cell of the cell list
for
(;
first_neigh_quad
!=
last_neigh_quad
;
++
first_neigh_quad
){
QuadraturePoint
quad
=
*
first_neigh_quad
;
UInt
nb_quad_per_elem
=
this
->
model
->
getFEEngine
().
getNbQuadraturePoints
(
quad
.
type
,
quad
.
ghost_type
);
// little optimization to not search in the map at each quad points
if
(
quad
.
type
!=
current_element_type
||
quad
.
ghost_type
!=
current_ghost_type
)
{
current_element_type
=
quad
.
type
;
current_ghost_type
=
quad
.
ghost_type
;
existing_pairs_num
=
quad
.
ghost_type
==
_not_ghost
?
0
:
1
;
if
(
!
pairs
.
exists
(
current_element_type
,
current_ghost_type
))
{
neighbors
=
&
(
pairs
.
alloc
(
0
,
2
,
current_element_type
,
current_ghost_type
));
}
else
{
neighbors
=
&
(
pairs
(
current_element_type
,
current_ghost_type
));
}
existing_pairs
[
existing_pairs_num
].
insert
(
std
::
pair
<
ElementType
,
ElementType
>
(
*
it
,
current_element_type
));
element_index_material_it
=
this
->
model
->
getElementIndexByMaterial
(
current_element_type
,
current_ghost_type
).
begin
(
2
);
quad_coord_it
=
quadrature_points_coordinates
(
current_element_type
,
current_ghost_type
).
begin
(
spatial_dimension
);
}
const
Vector
<
UInt
>
&
el_mat
=
element_index_material_it
[
quad
.
element
];
UInt
neigh_num_quad
=
el_mat
(
1
)
*
nb_quad_per_elem
+
quad
.
num_point
;
const
Vector
<
Real
>
&
neigh_quad
=
quad_coord_it
[
neigh_num_quad
];
Real
distance
=
first_quad
->
distance
(
neigh_quad
);
if
(
distance
<=
weight_func
->
getRadius
()
&&
(
quad
.
ghost_type
==
_ghost
||
(
quad
.
ghost_type
==
_not_ghost
&&
my_num_quad
<=
neigh_num_quad
)))
{
// storing only half lists
UInt
pair
[
2
];
pair
[
0
]
=
my_num_quad
;
pair
[
1
]
=
neigh_num_quad
;
neighbors
->
push_back
(
pair
);
// >>>>>> DEBUG CODE >>>>>> //
#if defined(AKANTU_DEBUG_TOOLS) && defined(AKANTU_CORE_CXX11)
debug
::
element_manager
.
print
(
debug
::
_dm_material_non_local
,
[
this
,
&
first_quad
,
&
neigh_quad
,
&
distance
,
&
quad_point
,
&
quad
](
const
Element
&
el
)
->
std
::
string
{
std
::
stringstream
out
;
if
((
Element
)
quad_point
==
el
)
{
out
<<
" neigh1: "
<<
quad
<<
" -- "
<<
*
first_quad
<<
" "
<<
neigh_quad
<<
" dist: "
<<
distance
;
}
if
((
Element
)
quad
==
el
)
{
out
<<
" neigh2: "
<<
quad_point
<<
" -- "
<<
*
first_quad
<<
" "
<<
neigh_quad
<<
" dist: "
<<
distance
;
}
return
out
.
str
();
});
#endif
// <<<<<< DEBUG CODE <<<<<< //
}
}
}
}
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
UInt
spatial_dimension
,
template
<
UInt
>
class
WeightFunction
>
void
MaterialNonLocal
<
spatial_dimension
,
WeightFunction
>::
computeWeights
(
const
ElementTypeMapArray
<
Real
>
&
quadrature_points_coordinates
)
{
AKANTU_DEBUG_IN
();
GhostType
ghost_type1
;
ghost_type1
=
_not_ghost
;
InternalField
<
Real
>
quadrature_points_volumes
(
"quadrature_points_volumes"
,
*
this
);
quadrature_points_volumes
.
initialize
(
1
);
const
FEEngine
&
fem
=
this
->
model
->
getFEEngine
();
weight_func
->
updateInternals
(
quadrature_points_volumes
);
for
(
UInt
gt
=
_not_ghost
;
gt
<=
_ghost
;
++
gt
)
{
GhostType
ghost_type2
=
(
GhostType
)
gt
;
UInt
existing_pairs_num
=
gt
-
_not_ghost
;
pair_type
::
iterator
first_pair_types
=
existing_pairs
[
existing_pairs_num
].
begin
();
pair_type
::
iterator
last_pair_types
=
existing_pairs
[
existing_pairs_num
].
end
();
// Compute the weights
for
(;
first_pair_types
!=
last_pair_types
;
++
first_pair_types
)
{
ElementType
type1
=
first_pair_types
->
first
;
ElementType
type2
=
first_pair_types
->
second
;
const
Array
<
UInt
>
&
pairs
=
pair_list
(
type1
,
ghost_type1
)(
type2
,
ghost_type2
);
std
::
string
ghost_id
=
""
;
if
(
ghost_type1
==
_ghost
)
ghost_id
=
":ghost"
;
ElementTypeMapArray
<
Real
>
&
weights_type_1
=
pair_weight
(
type1
,
ghost_type1
);
std
::
stringstream
sstr
;
sstr
<<
getID
()
<<
":pair_weight:"
<<
type1
<<
ghost_id
;
weights_type_1
.
setID
(
sstr
.
str
());
Array
<
Real
>
*
tmp_weight
=
NULL
;
if
(
!
weights_type_1
.
exists
(
type2
,
ghost_type2
))
{
tmp_weight
=
&
(
weights_type_1
.
alloc
(
0
,
2
,
type2
,
ghost_type2
));
}
else
{
tmp_weight
=
&
(
weights_type_1
(
type2
,
ghost_type2
));
}
Array
<
Real
>
&
weights
=
*
tmp_weight
;
weights
.
resize
(
pairs
.
getSize
());
weights
.
clear
();
const
Array
<
Real
>
&
jacobians_1
=
fem
.
getIntegratorInterface
().
getJacobians
(
type1
,
ghost_type1
);
const
Array
<
Real
>
&
jacobians_2
=
fem
.
getIntegratorInterface
().
getJacobians
(
type2
,
ghost_type2
);
const
Array
<
UInt
>
&
elem_mat_1
=
element_filter
(
type1
,
ghost_type1
);
const
Array
<
UInt
>
&
elem_mat_2
=
element_filter
(
type2
,
ghost_type2
);
UInt
nb_quad1
=
fem
.
getNbQuadraturePoints
(
type1
);
UInt
nb_quad2
=
fem
.
getNbQuadraturePoints
(
type2
);
Array
<
Real
>
&
quads_volumes1
=
quadrature_points_volumes
(
type1
,
ghost_type1
);
Array
<
Real
>
&
quads_volumes2
=
quadrature_points_volumes
(
type2
,
ghost_type2
);
Array
<
Real
>::
const_vector_iterator
iquads1
;
Array
<
Real
>::
const_vector_iterator
iquads2
;
iquads1
=
quadrature_points_coordinates
(
type1
,
ghost_type1
).
begin
(
spatial_dimension
);
iquads2
=
quadrature_points_coordinates
(
type2
,
ghost_type2
).
begin
(
spatial_dimension
);
Array
<
UInt
>::
const_iterator
<
Vector
<
UInt
>
>
first_pair
=
pairs
.
begin
(
2
);
Array
<
UInt
>::
const_iterator
<
Vector
<
UInt
>
>
last_pair
=
pairs
.
end
(
2
);
Array
<
Real
>::
vector_iterator
weight
=
weights
.
begin
(
2
);
this
->
weight_func
->
selectType
(
type1
,
ghost_type1
,
type2
,
ghost_type2
);
// Weight function
for
(;
first_pair
!=
last_pair
;
++
first_pair
,
++
weight
)
{
UInt
_q1
=
(
*
first_pair
)(
0
);
UInt
_q2
=
(
*
first_pair
)(
1
);
const
Vector
<
Real
>
&
pos1
=
iquads1
[
_q1
];
const
Vector
<
Real
>
&
pos2
=
iquads2
[
_q2
];
QuadraturePoint
q1
(
_q1
/
nb_quad1
,
_q1
%
nb_quad1
,
_q1
,
pos1
,
type1
,
ghost_type1
);
QuadraturePoint
q2
(
_q2
/
nb_quad2
,
_q2
%
nb_quad2
,
_q2
,
pos2
,
type2
,
ghost_type2
);
Real
r
=
pos1
.
distance
(
pos2
);
Real
w2J2
=
jacobians_2
(
elem_mat_2
(
q2
.
element
)
*
nb_quad2
+
q2
.
num_point
);
Real
w
=
this
->
weight_func
->
operator
()(
r
,
q1
,
q2
);
(
*
weight
)(
0
)
=
w2J2
*
w
;
quads_volumes1
(
_q1
)
+=
(
*
weight
)(
0
);
if
(
ghost_type2
!=
_ghost
&&
_q1
!=
_q2
)
{
Real
w1J1
=
jacobians_1
(
elem_mat_1
(
q1
.
element
)
*
nb_quad1
+
q1
.
num_point
);
(
*
weight
)(
1
)
=
w1J1
*
this
->
weight_func
->
operator
()(
r
,
q2
,
q1
);
quads_volumes2
(
_q2
)
+=
(
*
weight
)(
1
);
}
else
(
*
weight
)(
1
)
=
0
;
}
}
}
//normalize the weights
for
(
UInt
gt
=
_not_ghost
;
gt
<=
_ghost
;
++
gt
)
{
GhostType
ghost_type2
=
(
GhostType
)
gt
;
UInt
existing_pairs_num
=
gt
-
_not_ghost
;
pair_type
::
iterator
first_pair_types
=
existing_pairs
[
existing_pairs_num
].
begin
();
pair_type
::
iterator
last_pair_types
=
existing_pairs
[
existing_pairs_num
].
end
();
first_pair_types
=
existing_pairs
[
existing_pairs_num
].
begin
();
for
(;
first_pair_types
!=
last_pair_types
;
++
first_pair_types
)
{
ElementType
type1
=
first_pair_types
->
first
;
ElementType
type2
=
first_pair_types
->
second
;
const
Array
<
UInt
>
&
pairs
=
pair_list
(
type1
,
ghost_type1
)(
type2
,
ghost_type2
);
Array
<
Real
>
&
weights
=
pair_weight
(
type1
,
ghost_type1
)(
type2
,
ghost_type2
);
Array
<
Real
>
&
quads_volumes1
=
quadrature_points_volumes
(
type1
,
ghost_type1
);
Array
<
Real
>
&
quads_volumes2
=
quadrature_points_volumes
(
type2
,
ghost_type2
);
Array
<
UInt
>::
const_iterator
<
Vector
<
UInt
>
>
first_pair
=
pairs
.
begin
(
2
);
Array
<
UInt
>::
const_iterator
<
Vector
<
UInt
>
>
last_pair
=
pairs
.
end
(
2
);
Array
<
Real
>::
vector_iterator
weight
=
weights
.
begin
(
2
);
for
(;
first_pair
!=
last_pair
;
++
first_pair
,
++
weight
)
{
UInt
q1
=
(
*
first_pair
)(
0
);
UInt
q2
=
(
*
first_pair
)(
1
);
(
*
weight
)(
0
)
*=
1.
/
quads_volumes1
(
q1
);
if
(
ghost_type2
!=
_ghost
)
(
*
weight
)(
1
)
*=
1.
/
quads_volumes2
(
q2
);
}
}
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
UInt
spatial_dimension
,
template
<
UInt
>
class
WeightFunction
>
template
<
typename
T
>
void
MaterialNonLocal
<
spatial_dimension
,
WeightFunction
>::
weightedAvergageOnNeighbours
(
const
ElementTypeMapArray
<
T
>
&
to_accumulate
,
ElementTypeMapArray
<
T
>
&
accumulated
,
UInt
nb_degree_of_freedom
,
GhostType
ghost_type2
)
const
{
AKANTU_DEBUG_IN
();
UInt
existing_pairs_num
=
0
;
if
(
ghost_type2
==
_ghost
)
existing_pairs_num
=
1
;
pair_type
::
const_iterator
first_pair_types
=
existing_pairs
[
existing_pairs_num
].
begin
();
pair_type
::
const_iterator
last_pair_types
=
existing_pairs
[
existing_pairs_num
].
end
();
GhostType
ghost_type1
=
_not_ghost
;
// does not make sens the ghost vs ghost so this should always by _not_ghost
for
(;
first_pair_types
!=
last_pair_types
;
++
first_pair_types
)
{
const
Array
<
UInt
>
&
pairs
=
pair_list
(
first_pair_types
->
first
,
ghost_type1
)(
first_pair_types
->
second
,
ghost_type2
);
const
Array
<
Real
>
&
weights
=
pair_weight
(
first_pair_types
->
first
,
ghost_type1
)(
first_pair_types
->
second
,
ghost_type2
);
const
Array
<
T
>
&
to_acc
=
to_accumulate
(
first_pair_types
->
second
,
ghost_type2
);
Array
<
T
>
&
acc
=
accumulated
(
first_pair_types
->
first
,
ghost_type1
);
if
(
ghost_type2
==
_not_ghost
)
acc
.
clear
();
Array
<
UInt
>::
const_iterator
<
Vector
<
UInt
>
>
first_pair
=
pairs
.
begin
(
2
);
Array
<
UInt
>::
const_iterator
<
Vector
<
UInt
>
>
last_pair
=
pairs
.
end
(
2
);
Array
<
Real
>::
const_vector_iterator
pair_w
=
weights
.
begin
(
2
);
for
(;
first_pair
!=
last_pair
;
++
first_pair
,
++
pair_w
)
{
UInt
q1
=
(
*
first_pair
)(
0
);
UInt
q2
=
(
*
first_pair
)(
1
);
for
(
UInt
d
=
0
;
d
<
nb_degree_of_freedom
;
++
d
){
acc
(
q1
,
d
)
+=
(
*
pair_w
)(
0
)
*
to_acc
(
q2
,
d
);
if
(
ghost_type2
!=
_ghost
)
acc
(
q2
,
d
)
+=
(
*
pair_w
)(
1
)
*
to_acc
(
q1
,
d
);
}
}
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
UInt
spatial_dimension
,
template
<
UInt
>
class
WeightFunction
>
void
MaterialNonLocal
<
spatial_dimension
,
WeightFunction
>::
updateResidual
(
GhostType
ghost_type
)
{
AKANTU_DEBUG_IN
();
// Update the weights for the non local variable averaging
if
(
ghost_type
==
_not_ghost
&&
this
->
weight_func
->
getUpdateRate
()
&&
(
this
->
compute_stress_calls
%
this
->
weight_func
->
getUpdateRate
()
==
0
))
{
ElementTypeMapArray
<
Real
>
quadrature_points_coordinates
(
"quadrature_points_coordinates"
,
getID
());
Mesh
&
mesh
=
this
->
model
->
getFEEngine
().
getMesh
();
mesh
.
initElementTypeMapArray
(
quadrature_points_coordinates
,
spatial_dimension
,
spatial_dimension
);
computeQuadraturePointsCoordinates
(
quadrature_points_coordinates
,
_not_ghost
);
computeQuadraturePointsCoordinates
(
quadrature_points_coordinates
,
_ghost
);
computeWeights
(
quadrature_points_coordinates
);
}
if
(
ghost_type
==
_not_ghost
)
++
this
->
compute_stress_calls
;
computeAllStresses
(
ghost_type
);
computeNonLocalStresses
(
ghost_type
);
assembleResidual
(
ghost_type
);
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
UInt
spatial_dimension
,
template
<
UInt
>
class
WeightFunction
>
void
MaterialNonLocal
<
spatial_dimension
,
WeightFunction
>::
computeAllNonLocalStresses
(
GhostType
ghost_type
)
{
// Update the weights for the non local variable averaging
if
(
ghost_type
==
_not_ghost
)
{
if
(
this
->
weight_func
->
getUpdateRate
()
&&
(
this
->
compute_stress_calls
%
this
->
weight_func
->
getUpdateRate
()
==
0
))
{
this
->
model
->
getSynchronizerRegistry
().
asynchronousSynchronize
(
_gst_mnl_weight
);
ElementTypeMapArray
<
Real
>
quadrature_points_coordinates
(
"quadrature_points_coordinates"
,
getID
());
Mesh
&
mesh
=
this
->
model
->
getFEEngine
().
getMesh
();
mesh
.
initElementTypeMapArray
(
quadrature_points_coordinates
,
spatial_dimension
,
spatial_dimension
);
computeQuadraturePointsCoordinates
(
quadrature_points_coordinates
,
_not_ghost
);
computeQuadraturePointsCoordinates
(
quadrature_points_coordinates
,
_ghost
);
this
->
model
->
getSynchronizerRegistry
().
waitEndSynchronize
(
_gst_mnl_weight
);
computeWeights
(
quadrature_points_coordinates
);
}
typename
std
::
map
<
ID
,
NonLocalVariable
>::
iterator
it
=
non_local_variables
.
begin
();
typename
std
::
map
<
ID
,
NonLocalVariable
>::
iterator
end
=
non_local_variables
.
end
();
for
(;
it
!=
end
;
++
it
)
{
NonLocalVariable
&
non_local_variable
=
it
->
second
;
non_local_variable
.
non_local
->
resize
();
this
->
weightedAvergageOnNeighbours
(
*
non_local_variable
.
local
,
*
non_local_variable
.
non_local
,
non_local_variable
.
nb_component
,
_not_ghost
);
}
++
this
->
compute_stress_calls
;
}
else
{
typename
std
::
map
<
ID
,
NonLocalVariable
>::
iterator
it
=
non_local_variables
.
begin
();
typename
std
::
map
<
ID
,
NonLocalVariable
>::
iterator
end
=
non_local_variables
.
end
();
for
(;
it
!=
end
;
++
it
)
{
NonLocalVariable
&
non_local_variable
=
it
->
second
;
this
->
weightedAvergageOnNeighbours
(
*
non_local_variable
.
local
,
*
non_local_variable
.
non_local
,
non_local_variable
.
nb_component
,
_ghost
);
// >>>>>> DEBUG CODE >>>>>> //
#if defined(AKANTU_DEBUG_TOOLS)
#if defined(AKANTU_CORE_CXX11)
debug
::
element_manager
.
print
(
debug
::
_dm_material
,
[
non_local_variable
,
&
element_filter
](
const
Element
&
el
)
->
std
::
string
{
std
::
stringstream
out
;
if
(
element_filter
.
exists
(
el
.
type
,
el
.
ghost_type
))
{
UInt
pos
=
element_filter
(
el
.
type
,
el
.
ghost_type
).
find
(
el
.
element
);
if
(
pos
!=
UInt
(
-
1
))
{
out
<<
(
*
non_local_variable
.
local
)(
el
.
type
,
el
.
ghost_type
).
getID
()
<<
": loc "
<<
(
*
non_local_variable
.
local
)(
el
.
type
,
el
.
ghost_type
).
begin
(
non_local_variable
.
nb_component
)[
pos
]
<<
" - non loc "
<<
(
*
non_local_variable
.
non_local
)(
el
.
type
,
el
.
ghost_type
).
begin
(
non_local_variable
.
nb_component
)[
pos
];
}
}
return
out
.
str
();
});
#else
debug
::
element_manager
.
printData
(
debug
::
_dm_material
,
"MaterialNonLocal: computeAllNonLocalStress"
,
*
non_local_variable
.
local
,
this
->
element_filter
);
debug
::
element_manager
.
printData
(
debug
::
_dm_material
,
"MaterialNonLocal: computeAllNonLocalStress"
,
*
non_local_variable
.
non_local
,
this
->
element_filter
);
#endif
#endif
// <<<<<< DEBUG CODE <<<<<< //
}
computeNonLocalStresses
(
_not_ghost
);
// >>>>>> DEBUG CODE >>>>>> //
#if defined(AKANTU_DEBUG_TOOLS) && defined(AKANTU_CORE_CXX11)
ElementTypeMapArray
<
Real
>
quadrature_points_coordinates
(
"quadrature_points_coordinates"
,
id
);
Mesh
&
mesh
=
this
->
model
->
getFEEngine
().
getMesh
();
mesh
.
initElementTypeMapArray
(
quadrature_points_coordinates
,
spatial_dimension
,
spatial_dimension
);
computeQuadraturePointsCoordinates
(
quadrature_points_coordinates
,
_not_ghost
);
computeQuadraturePointsCoordinates
(
quadrature_points_coordinates
,
_ghost
);
debug
::
element_manager
.
print
(
debug
::
_dm_material_non_local
,
[
this
,
&
quadrature_points_coordinates
](
const
Element
&
el
)
->
std
::
string
{
static
UInt
step
=
0
;
std
::
stringstream
out
;
GhostType
ghost_type1
=
_not_ghost
;
FEEngine
&
fem
=
this
->
model
->
getFEEngine
();
Mesh
&
mesh
=
this
->
model
->
getMesh
();
std
::
ofstream
quad_out
;
std
::
stringstream
sstro
;
sstro
<<
"neigh_vals_"
<<
el
.
element
<<
".csv"
;
if
(
step
==
0
)
{
quad_out
.
open
(
sstro
.
str
());
quad_out
<<
"#id,step,gt,realid,w1,w2"
;
for
(
UInt
i
=
0
;
i
<
spatial_dimension
;
++
i
)
{
std
::
stringstream
sstr
;
sstr
<<
",x"
<<
i
;
quad_out
<<
sstr
.
str
();
}
quad_out
<<
",dam,jac"
;
typename
std
::
map
<
ID
,
NonLocalVariable
>::
iterator
it
=
non_local_variables
.
begin
();
typename
std
::
map
<
ID
,
NonLocalVariable
>::
iterator
end
=
non_local_variables
.
end
();
for
(;
it
!=
end
;
++
it
)
{
NonLocalVariable
&
non_local_variable
=
it
->
second
;
for
(
UInt
i
=
0
;
i
<
non_local_variable
.
nb_component
;
++
i
)
{
std
::
stringstream
sstr
;
sstr
<<
","
<<
non_local_variable
.
local
->
getID
()
<<
i
;
quad_out
<<
sstr
.
str
();
}
for
(
UInt
i
=
0
;
i
<
non_local_variable
.
nb_component
;
++
i
)
{
std
::
stringstream
sstr
;
sstr
<<
","
<<
non_local_variable
.
non_local
->
getID
()
<<
i
;
quad_out
<<
sstr
.
str
();
}
}
quad_out
<<
std
::
endl
;
}
else
{
quad_out
.
open
(
sstro
.
str
(),
std
::
ios_base
::
app
);
}
quad_out
.
precision
(
16
);
for
(
ghost_type_t
::
iterator
git
=
ghost_type_t
::
begin
();
git
!=
ghost_type_t
::
end
();
++
git
)
{
GhostType
ghost_type2
=
*
git
;
UInt
existing_pairs_num
=
ghost_type2
-
_not_ghost
;
pair_type
::
iterator
first_pair_types
=
existing_pairs
[
existing_pairs_num
].
begin
();
pair_type
::
iterator
last_pair_types
=
existing_pairs
[
existing_pairs_num
].
end
();
// Compute the weights
for
(;
first_pair_types
!=
last_pair_types
;
++
first_pair_types
)
{
ElementType
type1
=
first_pair_types
->
first
;
ElementType
type2
=
first_pair_types
->
second
;
const
Array
<
UInt
>
&
elem_mat_1
=
element_filter
(
type1
,
ghost_type1
);
const
Array
<
UInt
>
&
elem_mat_2
=
element_filter
(
type2
,
ghost_type2
);
UInt
nb_quad1
=
fem
.
getNbQuadraturePoints
(
type1
,
ghost_type1
);
UInt
nb_quad2
=
fem
.
getNbQuadraturePoints
(
type2
,
ghost_type2
);
const
Array
<
UInt
>
&
pairs
=
pair_list
(
type1
,
ghost_type1
)(
type2
,
ghost_type2
);
const
Array
<
Real
>
&
weights
=
pair_weight
(
type1
,
ghost_type1
)(
type2
,
ghost_type2
);
Array
<
UInt
>::
const_iterator
<
Vector
<
UInt
>
>
first_pair
=
pairs
.
begin
(
2
);
Array
<
UInt
>::
const_iterator
<
Vector
<
UInt
>
>
last_pair
=
pairs
.
end
(
2
);
Array
<
Real
>::
const_vector_iterator
pair_w
=
weights
.
begin
(
2
);
for
(;
first_pair
!=
last_pair
;
++
first_pair
,
++
pair_w
)
{
UInt
_q1
=
(
*
first_pair
)(
0
);
UInt
_q2
=
(
*
first_pair
)(
1
);
QuadraturePoint
q1
(
type1
,
elem_mat_1
(
_q1
/
nb_quad1
),
_q1
%
nb_quad1
,
ghost_type1
);
QuadraturePoint
q2
(
type2
,
elem_mat_2
(
_q2
/
nb_quad2
),
_q2
%
nb_quad2
,
ghost_type2
);
q1
.
global_num
=
_q1
;
q2
.
global_num
=
_q2
;
if
(
el
==
(
Element
)
q1
||
el
==
(
Element
)
q2
)
{
QuadraturePoint
q
;
Real
w1
,
w2
;
UInt
nb_quad
;
if
(
el
==
(
Element
)
q1
)
{
q
=
q2
;
nb_quad
=
nb_quad2
;
w1
=
(
*
pair_w
)(
0
);
w2
=
(
*
pair_w
)(
1
);
}
else
{
q
=
q1
;
nb_quad
=
nb_quad1
;
w1
=
(
*
pair_w
)(
1
);
w2
=
(
*
pair_w
)(
0
);
}
quad_out
<<
((
q
.
ghost_type
==
_ghost
?
mesh
.
getNbElement
(
q
.
type
,
q
.
ghost_type
)
:
0
)
+
q
.
global_num
)
<<
","
<<
step
<<
","
<<
q
.
ghost_type
<<
","
<<
q
.
global_num
;
quad_out
<<
","
<<
w1
<<
","
<<
w2
;
Array
<
Real
>::
const_vector_iterator
iquads
=
quadrature_points_coordinates
(
q
.
type
,
q
.
ghost_type
).
begin
(
spatial_dimension
);
const
Vector
<
Real
>
&
coord
=
iquads
[
q
.
global_num
];
for
(
UInt
i
(
0
);
i
<
coord
.
size
();
++
i
)
{
quad_out
<<
","
<<
coord
(
i
);
}
typename
std
::
map
<
ID
,
NonLocalVariable
>::
iterator
nlit
=
non_local_variables
.
begin
();
typename
std
::
map
<
ID
,
NonLocalVariable
>::
iterator
nlend
=
non_local_variables
.
end
();
for
(;
nlit
!=
nlend
;
++
nlit
)
{
NonLocalVariable
&
non_local_variable
=
nlit
->
second
;
const
Array
<
Real
>
&
local
=
(
*
non_local_variable
.
local
)(
q
.
type
,
q
.
ghost_type
);
for
(
UInt
i
(
0
);
i
<
local
.
getNbComponent
();
++
i
)
{
quad_out
<<
","
<<
local
(
q
.
global_num
,
i
);
}
const
Array
<
Real
>
&
non_local
=
(
*
non_local_variable
.
non_local
)(
q
.
type
,
q
.
ghost_type
);
for
(
UInt
i
(
0
);
i
<
non_local
.
getNbComponent
();
++
i
)
{
quad_out
<<
","
<<
non_local
(
q
.
global_num
,
i
);
}
}
Real
d
=
getArray
(
"damage"
,
q
.
type
,
q
.
ghost_type
)(
q
.
global_num
);
Real
j
=
fem
.
getIntegratorInterface
().
getJacobians
(
q
.
type
,
q
.
ghost_type
)(
q
.
element
*
nb_quad
+
q
.
num_point
);
quad_out
<<
","
<<
d
<<
","
<<
j
;
quad_out
<<
std
::
endl
;
}
}
}
}
step
++
;
return
out
.
str
();
});
#endif
// <<<<<< DEBUG CODE <<<<<< //
}
// >>>>>> DEBUG CODE >>>>>> //
#if defined(AKANTU_DEBUG_TOOLS)
#if defined(AKANTU_CORE_CXX11)
debug
::
element_manager
.
print
(
debug
::
_dm_material_damage
,
[
ghost_type
,
this
](
const
Element
&
el
)
->
std
::
string
{
std
::
stringstream
out
;
if
(
el
.
ghost_type
==
ghost_type
&&
element_filter
.
exists
(
el
.
type
,
ghost_type
))
{
UInt
pos
=
element_filter
(
el
.
type
,
el
.
ghost_type
).
find
(
el
.
element
);
if
(
pos
!=
UInt
(
-
1
))
{
Real
d
=
getArray
(
"damage"
,
el
.
type
,
el
.
ghost_type
)(
pos
);
out
<<
" damage: "
<<
d
;
}
}
return
out
.
str
();
});
#endif
#endif
// <<<<<< DEBUG CODE <<<<<< //
}
/* -------------------------------------------------------------------------- */
template
<
UInt
spatial_dimension
,
template
<
UInt
>
class
WeightFunction
>
void
MaterialNonLocal
<
spatial_dimension
,
WeightFunction
>::
savePairs
(
const
std
::
string
&
filename
)
const
{
std
::
ofstream
pout
;
std
::
stringstream
sstr
;
StaticCommunicator
&
comm
=
StaticCommunicator
::
getStaticCommunicator
();
Int
prank
=
comm
.
whoAmI
();
sstr
<<
filename
<<
"."
<<
prank
;
pout
.
open
(
sstr
.
str
().
c_str
());
GhostType
ghost_type1
;
ghost_type1
=
_not_ghost
;
for
(
UInt
gt
=
_not_ghost
;
gt
<=
_ghost
;
++
gt
)
{
GhostType
ghost_type2
=
(
GhostType
)
gt
;
UInt
existing_pairs_num
=
gt
-
_not_ghost
;
pair_type
::
const_iterator
first_pair_types
=
existing_pairs
[
existing_pairs_num
].
begin
();
pair_type
::
const_iterator
last_pair_types
=
existing_pairs
[
existing_pairs_num
].
end
();
for
(;
first_pair_types
!=
last_pair_types
;
++
first_pair_types
)
{
const
Array
<
UInt
>
&
pairs
=
(
*
pair_list
(
first_pair_types
->
first
,
ghost_type1
))(
first_pair_types
->
second
,
ghost_type2
);
const
Array
<
Real
>
&
weights
=
(
*
pair_weight
(
first_pair_types
->
first
,
ghost_type1
))(
first_pair_types
->
second
,
ghost_type2
);
pout
<<
"Types : "
<<
first_pair_types
->
first
<<
" ("
<<
ghost_type1
<<
") - "
<<
first_pair_types
->
second
<<
" ("
<<
ghost_type2
<<
")"
<<
std
::
endl
;
Array
<
UInt
>::
const_iterator
<
Vector
<
UInt
>
>
first_pair
=
pairs
.
begin
(
2
);
Array
<
UInt
>::
const_iterator
<
Vector
<
UInt
>
>
last_pair
=
pairs
.
end
(
2
);
Array
<
Real
>::
const_vector_iterator
pair_w
=
weights
.
begin
(
2
);
for
(;
first_pair
!=
last_pair
;
++
first_pair
,
++
pair_w
)
{
UInt
q1
=
(
*
first_pair
)(
0
);
UInt
q2
=
(
*
first_pair
)(
1
);
pout
<<
q1
<<
" "
<<
q2
<<
" "
<<
(
*
pair_w
)(
0
)
<<
" "
<<
(
*
pair_w
)(
1
)
<<
std
::
endl
;
}
}
}
}
/* -------------------------------------------------------------------------- */
template
<
UInt
spatial_dimension
,
template
<
UInt
>
class
WeightFunction
>
void
MaterialNonLocal
<
spatial_dimension
,
WeightFunction
>::
neighbourhoodStatistics
(
const
std
::
string
&
filename
)
const
{
std
::
ofstream
pout
;
pout
.
open
(
filename
.
c_str
());
const
Mesh
&
mesh
=
this
->
model
->
getFEEngine
().
getMesh
();
GhostType
ghost_type1
;
ghost_type1
=
_not_ghost
;
StaticCommunicator
&
comm
=
StaticCommunicator
::
getStaticCommunicator
();
Int
prank
=
comm
.
whoAmI
();
InternalField
<
UInt
>
nb_neighbors
(
"nb_neighbours"
,
*
const_cast
<
MaterialNonLocal
*>
(
this
));
nb_neighbors
.
initialize
(
1
);
for
(
UInt
gt
=
_not_ghost
;
gt
<=
_ghost
;
++
gt
)
{
GhostType
ghost_type2
=
(
GhostType
)
gt
;
UInt
existing_pairs_num
=
gt
-
_not_ghost
;
pair_type
::
const_iterator
first_pair_types
=
existing_pairs
[
existing_pairs_num
].
begin
();
pair_type
::
const_iterator
last_pair_types
=
existing_pairs
[
existing_pairs_num
].
end
();
for
(;
first_pair_types
!=
last_pair_types
;
++
first_pair_types
)
{
const
Array
<
UInt
>
&
pairs
=
pair_list
(
first_pair_types
->
first
,
ghost_type1
)(
first_pair_types
->
second
,
ghost_type2
);
if
(
prank
==
0
)
{
pout
<<
ghost_type2
<<
" "
;
pout
<<
"Types : "
<<
first_pair_types
->
first
<<
" "
<<
first_pair_types
->
second
<<
std
::
endl
;
}
Array
<
UInt
>::
const_iterator
<
Vector
<
UInt
>
>
first_pair
=
pairs
.
begin
(
2
);
Array
<
UInt
>::
const_iterator
<
Vector
<
UInt
>
>
last_pair
=
pairs
.
end
(
2
);
Array
<
UInt
>
&
nb_neigh_1
=
nb_neighbors
(
first_pair_types
->
first
,
ghost_type1
);
Array
<
UInt
>
&
nb_neigh_2
=
nb_neighbors
(
first_pair_types
->
second
,
ghost_type2
);
for
(;
first_pair
!=
last_pair
;
++
first_pair
)
{
UInt
q1
=
(
*
first_pair
)(
0
);
UInt
q2
=
(
*
first_pair
)(
1
);
++
(
nb_neigh_1
(
q1
));
if
(
q1
!=
q2
)
++
(
nb_neigh_2
(
q2
));
}
}
Mesh
::
type_iterator
it
=
mesh
.
firstType
(
spatial_dimension
,
ghost_type1
);
Mesh
::
type_iterator
last_type
=
mesh
.
lastType
(
spatial_dimension
,
ghost_type1
);
UInt
nb_quads
=
0
;
Real
sum_nb_neig
=
0
;
UInt
max_nb_neig
=
0
;
UInt
min_nb_neig
=
std
::
numeric_limits
<
UInt
>::
max
();
for
(;
it
!=
last_type
;
++
it
)
{
Array
<
UInt
>
&
nb_neighor
=
nb_neighbors
(
*
it
,
ghost_type1
);
Array
<
UInt
>::
iterator
<
UInt
>
nb_neigh
=
nb_neighor
.
begin
();
Array
<
UInt
>::
iterator
<
UInt
>
end_neigh
=
nb_neighor
.
end
();
for
(;
nb_neigh
!=
end_neigh
;
++
nb_neigh
,
++
nb_quads
)
{
UInt
nb
=
*
nb_neigh
;
sum_nb_neig
+=
nb
;
max_nb_neig
=
std
::
max
(
max_nb_neig
,
nb
);
min_nb_neig
=
std
::
min
(
min_nb_neig
,
nb
);
}
}
comm
.
allReduce
(
&
nb_quads
,
1
,
_so_sum
);
comm
.
allReduce
(
&
sum_nb_neig
,
1
,
_so_sum
);
comm
.
allReduce
(
&
max_nb_neig
,
1
,
_so_max
);
comm
.
allReduce
(
&
min_nb_neig
,
1
,
_so_min
);
if
(
prank
==
0
)
{
pout
<<
ghost_type2
<<
" "
;
pout
<<
"Nb quadrature points: "
<<
nb_quads
<<
std
::
endl
;
Real
mean_nb_neig
=
sum_nb_neig
/
Real
(
nb_quads
);
pout
<<
ghost_type2
<<
" "
;
pout
<<
"Average nb neighbors: "
<<
mean_nb_neig
<<
"("
<<
sum_nb_neig
<<
")"
<<
std
::
endl
;
pout
<<
ghost_type2
<<
" "
;
pout
<<
"Max nb neighbors: "
<<
max_nb_neig
<<
std
::
endl
;
pout
<<
ghost_type2
<<
" "
;
pout
<<
"Min nb neighbors: "
<<
min_nb_neig
<<
std
::
endl
;
}
}
pout
.
close
();
}
/* -------------------------------------------------------------------------- */
template
<
UInt
spatial_dimension
,
template
<
UInt
>
class
WeightFunction
>
inline
UInt
MaterialNonLocal
<
spatial_dimension
,
WeightFunction
>::
getNbDataForElements
(
const
Array
<
Element
>
&
elements
,
SynchronizationTag
tag
)
const
{
UInt
nb_quadrature_points
=
this
->
getModel
().
getNbQuadraturePoints
(
elements
);
UInt
size
=
0
;
if
(
tag
==
_gst_mnl_for_average
)
{
typename
std
::
map
<
ID
,
NonLocalVariable
>::
const_iterator
it
=
non_local_variables
.
begin
();
typename
std
::
map
<
ID
,
NonLocalVariable
>::
const_iterator
end
=
non_local_variables
.
end
();
for
(;
it
!=
end
;
++
it
)
{
const
NonLocalVariable
&
non_local_variable
=
it
->
second
;
size
+=
non_local_variable
.
nb_component
*
sizeof
(
Real
)
*
nb_quadrature_points
;
}
}
size
+=
weight_func
->
getNbDataForElements
(
elements
,
tag
);
return
size
;
}
/* -------------------------------------------------------------------------- */
template
<
UInt
spatial_dimension
,
template
<
UInt
>
class
WeightFunction
>
inline
void
MaterialNonLocal
<
spatial_dimension
,
WeightFunction
>::
packElementData
(
CommunicationBuffer
&
buffer
,
const
Array
<
Element
>
&
elements
,
SynchronizationTag
tag
)
const
{
if
(
tag
==
_gst_mnl_for_average
)
{
typename
std
::
map
<
ID
,
NonLocalVariable
>::
const_iterator
it
=
non_local_variables
.
begin
();
typename
std
::
map
<
ID
,
NonLocalVariable
>::
const_iterator
end
=
non_local_variables
.
end
();
for
(;
it
!=
end
;
++
it
)
{
const
NonLocalVariable
&
non_local_variable
=
it
->
second
;
this
->
packElementDataHelper
(
*
non_local_variable
.
local
,
buffer
,
elements
);
}
}
weight_func
->
packElementData
(
buffer
,
elements
,
tag
);
}
/* -------------------------------------------------------------------------- */
template
<
UInt
spatial_dimension
,
template
<
UInt
>
class
WeightFunction
>
inline
void
MaterialNonLocal
<
spatial_dimension
,
WeightFunction
>::
unpackElementData
(
CommunicationBuffer
&
buffer
,
const
Array
<
Element
>
&
elements
,
SynchronizationTag
tag
)
{
if
(
tag
==
_gst_mnl_for_average
)
{
typename
std
::
map
<
ID
,
NonLocalVariable
>::
iterator
it
=
non_local_variables
.
begin
();
typename
std
::
map
<
ID
,
NonLocalVariable
>::
iterator
end
=
non_local_variables
.
end
();
for
(;
it
!=
end
;
++
it
)
{
NonLocalVariable
&
non_local_variable
=
it
->
second
;
this
->
unpackElementDataHelper
(
*
non_local_variable
.
local
,
buffer
,
elements
);
}
}
weight_func
->
unpackElementData
(
buffer
,
elements
,
tag
);
}
/* -------------------------------------------------------------------------- */
// template<UInt spatial_dimension, template <UInt> class WeightFunction>
// inline void MaterialNonLocal<spatial_dimension, WeightFunction>::onElementsAdded(const Array<Element> & element_list) {
// AKANTU_DEBUG_IN();
// Material::onElementsAdded(element_list, event);
// AKANTU_DEBUG_OUT();
// }
/* -------------------------------------------------------------------------- */
template
<
UInt
spatial_dimension
,
template
<
UInt
>
class
WeightFunction
>
inline
void
MaterialNonLocal
<
spatial_dimension
,
WeightFunction
>::
onElementsRemoved
(
const
Array
<
Element
>
&
element_list
,
const
ElementTypeMapArray
<
UInt
>
&
new_numbering
,
__attribute__
((
unused
))
const
RemovedElementsEvent
&
event
)
{
AKANTU_DEBUG_IN
();
Material
::
onElementsRemoved
(
element_list
,
new_numbering
,
event
);
pair_type
::
const_iterator
first_pair_types
=
existing_pairs
[
1
].
begin
();
pair_type
::
const_iterator
last_pair_types
=
existing_pairs
[
1
].
end
();
// Renumber element to keep
for
(;
first_pair_types
!=
last_pair_types
;
++
first_pair_types
)
{
ElementType
type2
=
first_pair_types
->
second
;
GhostType
ghost_type2
=
_ghost
;
UInt
nb_quad2
=
this
->
model
->
getFEEngine
().
getNbQuadraturePoints
(
type2
);
Array
<
UInt
>
&
pairs
=
pair_list
(
first_pair_types
->
first
,
_not_ghost
)(
first_pair_types
->
second
,
ghost_type2
);
Array
<
UInt
>::
iterator
<
Vector
<
UInt
>
>
first_pair
=
pairs
.
begin
(
2
);
Array
<
UInt
>::
iterator
<
Vector
<
UInt
>
>
last_pair
=
pairs
.
end
(
2
);
for
(;
first_pair
!=
last_pair
;
++
first_pair
)
{
UInt
_q2
=
(
*
first_pair
)(
1
);
const
Array
<
UInt
>
&
renumbering
=
new_numbering
(
type2
,
ghost_type2
);
UInt
el
=
_q2
/
nb_quad2
;
UInt
new_el
=
renumbering
(
el
);
AKANTU_DEBUG_ASSERT
(
new_el
!=
UInt
(
-
1
),
"A local quad as been removed instead f just renumbered"
);
(
*
first_pair
)(
1
)
=
new_el
*
nb_quad2
+
_q2
%
nb_quad2
;
}
}
AKANTU_DEBUG_OUT
();
}
Event Timeline
Log In to Comment