Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F79599791
solid_mechanics_model_RVE.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
Tue, Aug 27, 07:22
Size
21 KB
Mime Type
text/x-c
Expires
Thu, Aug 29, 07:22 (2 d)
Engine
blob
Format
Raw Data
Handle
20208830
Attached To
rAKA akantu
solid_mechanics_model_RVE.cc
View Options
/**
* Copyright (©) 2018-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
*
* 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 "solid_mechanics_model_RVE.hh"
#include "element_group.hh"
#include "material_damage_iterative.hh"
#include "node_group.hh"
#include "non_linear_solver.hh"
#include "non_local_manager.hh"
#include "parser.hh"
#include "sparse_matrix.hh"
/* -------------------------------------------------------------------------- */
#include <string>
/* -------------------------------------------------------------------------- */
namespace
akantu
{
/* -------------------------------------------------------------------------- */
SolidMechanicsModelRVE
::
SolidMechanicsModelRVE
(
Mesh
&
mesh
,
bool
use_RVE_mat_selector
,
UInt
nb_gel_pockets
,
Int
dim
,
const
ID
&
id
)
:
SolidMechanicsModel
(
mesh
,
dim
,
id
),
volume
(
0.
),
use_RVE_mat_selector
(
use_RVE_mat_selector
),
nb_gel_pockets
(
nb_gel_pockets
),
nb_dumps
(
0
)
{
AKANTU_DEBUG_IN
();
/// find the four corner nodes of the RVE
findCornerNodes
();
/// remove the corner nodes from the surface node groups:
/// This most be done because corner nodes a not periodic
mesh
.
getElementGroup
(
"top"
).
removeNode
(
corner_nodes
(
2
));
mesh
.
getElementGroup
(
"top"
).
removeNode
(
corner_nodes
(
3
));
mesh
.
getElementGroup
(
"left"
).
removeNode
(
corner_nodes
(
3
));
mesh
.
getElementGroup
(
"left"
).
removeNode
(
corner_nodes
(
0
));
mesh
.
getElementGroup
(
"bottom"
).
removeNode
(
corner_nodes
(
1
));
mesh
.
getElementGroup
(
"bottom"
).
removeNode
(
corner_nodes
(
0
));
mesh
.
getElementGroup
(
"right"
).
removeNode
(
corner_nodes
(
2
));
mesh
.
getElementGroup
(
"right"
).
removeNode
(
corner_nodes
(
1
));
const
auto
&
bottom
=
mesh
.
getElementGroup
(
"bottom"
).
getNodeGroup
();
bottom_nodes
.
insert
(
bottom
.
begin
(),
bottom
.
end
());
const
auto
&
left
=
mesh
.
getElementGroup
(
"left"
).
getNodeGroup
();
left_nodes
.
insert
(
left
.
begin
(),
left
.
end
());
// /// enforce periodicity on the displacement fluctuations
// auto surface_pair_1 = std::make_pair("top", "bottom");
// auto surface_pair_2 = std::make_pair("right", "left");
// SurfacePairList surface_pairs_list;
// surface_pairs_list.push_back(surface_pair_1);
// surface_pairs_list.push_back(surface_pair_2);
// TODO: To Nicolas correct the PBCs
// this->setPBC(surface_pairs_list);
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
SolidMechanicsModelRVE
::~
SolidMechanicsModelRVE
()
=
default
;
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModelRVE
::
initFullImpl
(
const
ModelOptions
&
options
)
{
AKANTU_DEBUG_IN
();
auto
options_cp
(
options
);
options_cp
.
analysis_method
=
AnalysisMethod
::
_static
;
SolidMechanicsModel
::
initFullImpl
(
options_cp
);
// this->initMaterials();
auto
&
fem
=
this
->
getFEEngine
(
"SolidMechanicsFEEngine"
);
/// compute the volume of the RVE
GhostType
gt
=
_not_ghost
;
for
(
auto
element_type
:
this
->
mesh
.
elementTypes
(
spatial_dimension
,
gt
,
_ek_not_defined
))
{
Array
<
Real
>
Volume
(
this
->
mesh
.
getNbElement
(
element_type
)
*
fem
.
getNbIntegrationPoints
(
element_type
),
1
,
1.
);
this
->
volume
=
fem
.
integrate
(
Volume
,
element_type
);
}
std
::
cout
<<
"The volume of the RVE is "
<<
this
->
volume
<<
std
::
endl
;
/// dumping
std
::
stringstream
base_name
;
base_name
<<
this
->
id
;
this
->
setBaseName
(
base_name
.
str
());
this
->
addDumpFieldVector
(
"displacement"
);
this
->
addDumpField
(
"stress"
);
this
->
addDumpField
(
"grad_u"
);
this
->
addDumpField
(
"eigen_grad_u"
);
this
->
addDumpField
(
"blocked_dofs"
);
this
->
addDumpField
(
"material_index"
);
this
->
addDumpField
(
"damage"
);
this
->
addDumpField
(
"Sc"
);
this
->
addDumpField
(
"external_force"
);
this
->
addDumpField
(
"equivalent_stress"
);
this
->
addDumpField
(
"internal_force"
);
this
->
addDumpField
(
"delta_T"
);
this
->
dump
();
this
->
nb_dumps
+=
1
;
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModelRVE
::
applyBoundaryConditions
(
const
Matrix
<
Real
>
&
displacement_gradient
)
{
AKANTU_DEBUG_IN
();
/// get the position of the nodes
const
Array
<
Real
>
&
pos
=
mesh
.
getNodes
();
/// storage for the coordinates of a given node and the displacement that will
/// be applied
Vector
<
Real
>
x
(
spatial_dimension
);
Vector
<
Real
>
appl_disp
(
spatial_dimension
);
/// fix top right node
UInt
node
=
this
->
corner_nodes
(
2
);
x
(
0
)
=
pos
(
node
,
0
);
x
(
1
)
=
pos
(
node
,
1
);
appl_disp
.
mul
<
false
>
(
displacement_gradient
,
x
);
(
*
this
->
blocked_dofs
)(
node
,
0
)
=
true
;
(
*
this
->
displacement
)(
node
,
0
)
=
appl_disp
(
0
);
(
*
this
->
blocked_dofs
)(
node
,
1
)
=
true
;
(
*
this
->
displacement
)(
node
,
1
)
=
appl_disp
(
1
);
// (*this->blocked_dofs)(node,0) = true; (*this->displacement)(node,0) = 0.;
// (*this->blocked_dofs)(node,1) = true; (*this->displacement)(node,1) = 0.;
/// apply Hx at all the other corner nodes; H: displ. gradient
node
=
this
->
corner_nodes
(
0
);
x
(
0
)
=
pos
(
node
,
0
);
x
(
1
)
=
pos
(
node
,
1
);
appl_disp
.
mul
<
false
>
(
displacement_gradient
,
x
);
(
*
this
->
blocked_dofs
)(
node
,
0
)
=
true
;
(
*
this
->
displacement
)(
node
,
0
)
=
appl_disp
(
0
);
(
*
this
->
blocked_dofs
)(
node
,
1
)
=
true
;
(
*
this
->
displacement
)(
node
,
1
)
=
appl_disp
(
1
);
node
=
this
->
corner_nodes
(
1
);
x
(
0
)
=
pos
(
node
,
0
);
x
(
1
)
=
pos
(
node
,
1
);
appl_disp
.
mul
<
false
>
(
displacement_gradient
,
x
);
(
*
this
->
blocked_dofs
)(
node
,
0
)
=
true
;
(
*
this
->
displacement
)(
node
,
0
)
=
appl_disp
(
0
);
(
*
this
->
blocked_dofs
)(
node
,
1
)
=
true
;
(
*
this
->
displacement
)(
node
,
1
)
=
appl_disp
(
1
);
node
=
this
->
corner_nodes
(
3
);
x
(
0
)
=
pos
(
node
,
0
);
x
(
1
)
=
pos
(
node
,
1
);
appl_disp
.
mul
<
false
>
(
displacement_gradient
,
x
);
(
*
this
->
blocked_dofs
)(
node
,
0
)
=
true
;
(
*
this
->
displacement
)(
node
,
0
)
=
appl_disp
(
0
);
(
*
this
->
blocked_dofs
)(
node
,
1
)
=
true
;
(
*
this
->
displacement
)(
node
,
1
)
=
appl_disp
(
1
);
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModelRVE
::
applyHomogeneousTemperature
(
const
Real
&
temperature
)
{
for
(
UInt
m
=
0
;
m
<
this
->
getNbMaterials
();
++
m
)
{
Material
&
mat
=
this
->
getMaterial
(
m
);
const
ElementTypeMapArray
<
Idx
>
&
filter_map
=
mat
.
getElementFilter
();
// Loop over all element types
for
(
auto
&&
type
:
filter_map
.
elementTypes
(
spatial_dimension
))
{
const
Array
<
Idx
>
&
filter
=
filter_map
(
type
);
if
(
filter
.
size
()
==
0
)
continue
;
Array
<
Real
>
&
delta_T
=
mat
.
getArray
<
Real
>
(
"delta_T"
,
type
);
Array
<
Real
>::
scalar_iterator
delta_T_it
=
delta_T
.
begin
();
Array
<
Real
>::
scalar_iterator
delta_T_end
=
delta_T
.
end
();
for
(;
delta_T_it
!=
delta_T_end
;
++
delta_T_it
)
{
*
delta_T_it
=
temperature
;
}
}
}
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModelRVE
::
findCornerNodes
()
{
AKANTU_DEBUG_IN
();
// find corner nodes
const
auto
&
position
=
mesh
.
getNodes
();
const
auto
&
lower_bounds
=
mesh
.
getLowerBounds
();
const
auto
&
upper_bounds
=
mesh
.
getUpperBounds
();
AKANTU_DEBUG_ASSERT
(
spatial_dimension
==
2
,
"This is 2D only!"
);
corner_nodes
.
resize
(
4
);
corner_nodes
.
set
(
UInt
(
-
1
));
for
(
auto
&&
data
:
enumerate
(
make_view
(
position
,
spatial_dimension
)))
{
auto
node
=
std
::
get
<
0
>
(
data
);
const
auto
&
X
=
std
::
get
<
1
>
(
data
);
auto
distance
=
X
.
distance
(
lower_bounds
);
// node 1
if
(
Math
::
are_float_equal
(
distance
,
0
))
{
corner_nodes
(
0
)
=
node
;
}
// node 2
else
if
(
Math
::
are_float_equal
(
X
(
_x
),
upper_bounds
(
_x
))
&&
Math
::
are_float_equal
(
X
(
_y
),
lower_bounds
(
_y
)))
{
corner_nodes
(
1
)
=
node
;
}
// node 3
else
if
(
Math
::
are_float_equal
(
X
(
_x
),
upper_bounds
(
_x
))
&&
Math
::
are_float_equal
(
X
(
_y
),
upper_bounds
(
_y
)))
{
corner_nodes
(
2
)
=
node
;
}
// node 4
else
if
(
Math
::
are_float_equal
(
X
(
_x
),
lower_bounds
(
_x
))
&&
Math
::
are_float_equal
(
X
(
_y
),
upper_bounds
(
_y
)))
{
corner_nodes
(
3
)
=
node
;
}
}
for
(
UInt
i
=
0
;
i
<
corner_nodes
.
size
();
++
i
)
{
if
(
corner_nodes
(
i
)
==
UInt
(
-
1
))
AKANTU_ERROR
(
"The corner node "
<<
i
+
1
<<
" wasn't found"
);
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModelRVE
::
advanceASR
(
const
Matrix
<
Real
>
&
prestrain
)
{
AKANTU_DEBUG_IN
();
AKANTU_DEBUG_ASSERT
(
spatial_dimension
==
2
,
"This is 2D only!"
);
/// apply the new eigenstrain
for
(
auto
element_type
:
mesh
.
elementTypes
(
spatial_dimension
,
_not_ghost
,
_ek_not_defined
))
{
Array
<
Real
>
&
prestrain_vect
=
const_cast
<
Array
<
Real
>
&>
(
this
->
getMaterial
(
"gel"
).
getInternal
<
Real
>
(
"eigen_grad_u"
)(
element_type
));
auto
prestrain_it
=
prestrain_vect
.
begin
(
spatial_dimension
,
spatial_dimension
);
auto
prestrain_end
=
prestrain_vect
.
end
(
spatial_dimension
,
spatial_dimension
);
for
(;
prestrain_it
!=
prestrain_end
;
++
prestrain_it
)
(
*
prestrain_it
)
=
prestrain
;
}
/// advance the damage
MaterialDamageIterative
<
2
>
&
mat_paste
=
dynamic_cast
<
MaterialDamageIterative
<
2
>
&>
(
*
this
->
materials
[
1
]);
MaterialDamageIterative
<
2
>
&
mat_aggregate
=
dynamic_cast
<
MaterialDamageIterative
<
2
>
&>
(
*
this
->
materials
[
0
]);
UInt
nb_damaged_elements
=
0
;
Real
max_eq_stress_aggregate
=
0
;
Real
max_eq_stress_paste
=
0
;
auto
&
solver
=
this
->
getNonLinearSolver
();
solver
.
set
(
"max_iterations"
,
2
);
solver
.
set
(
"threshold"
,
1e-6
);
solver
.
set
(
"convergence_type"
,
SolveConvergenceCriteria
::
_solution
);
do
{
this
->
solveStep
();
/// compute damage
max_eq_stress_aggregate
=
mat_aggregate
.
getNormMaxEquivalentStress
();
max_eq_stress_paste
=
mat_paste
.
getNormMaxEquivalentStress
();
nb_damaged_elements
=
0
;
if
(
max_eq_stress_aggregate
>
max_eq_stress_paste
)
nb_damaged_elements
=
mat_aggregate
.
updateDamage
();
else
if
(
max_eq_stress_aggregate
<
max_eq_stress_paste
)
nb_damaged_elements
=
mat_paste
.
updateDamage
();
else
nb_damaged_elements
=
(
mat_paste
.
updateDamage
()
+
mat_aggregate
.
updateDamage
());
std
::
cout
<<
"the number of damaged elements is "
<<
nb_damaged_elements
<<
std
::
endl
;
}
while
(
nb_damaged_elements
);
if
(
this
->
nb_dumps
%
10
==
0
)
{
this
->
dump
();
}
this
->
nb_dumps
+=
1
;
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
Real
SolidMechanicsModelRVE
::
averageTensorField
(
UInt
row_index
,
UInt
col_index
,
const
ID
&
field_type
)
{
AKANTU_DEBUG_IN
();
auto
&
fem
=
this
->
getFEEngine
(
"SolidMechanicsFEEngine"
);
Real
average
=
0
;
GhostType
gt
=
_not_ghost
;
for
(
auto
element_type
:
mesh
.
elementTypes
(
spatial_dimension
,
gt
,
_ek_not_defined
))
{
if
(
field_type
==
"stress"
)
{
for
(
UInt
m
=
0
;
m
<
this
->
materials
.
size
();
++
m
)
{
const
auto
&
stress_vec
=
this
->
materials
[
m
]
->
getStress
(
element_type
);
const
auto
&
elem_filter
=
this
->
materials
[
m
]
->
getElementFilter
(
element_type
);
Array
<
Real
>
int_stress_vec
(
elem_filter
.
size
(),
spatial_dimension
*
spatial_dimension
,
"int_of_stress"
);
fem
.
integrate
(
stress_vec
,
int_stress_vec
,
spatial_dimension
*
spatial_dimension
,
element_type
,
_not_ghost
,
elem_filter
);
for
(
UInt
k
=
0
;
k
<
elem_filter
.
size
();
++
k
)
average
+=
int_stress_vec
(
k
,
row_index
*
spatial_dimension
+
col_index
);
// 3 is the value
// for the yy (in
// 3D, the value is
// 4)
}
}
else
if
(
field_type
==
"strain"
)
{
for
(
UInt
m
=
0
;
m
<
this
->
materials
.
size
();
++
m
)
{
const
auto
&
gradu_vec
=
this
->
materials
[
m
]
->
getGradU
(
element_type
);
const
auto
&
elem_filter
=
this
->
materials
[
m
]
->
getElementFilter
(
element_type
);
Array
<
Real
>
int_gradu_vec
(
elem_filter
.
size
(),
spatial_dimension
*
spatial_dimension
,
"int_of_gradu"
);
fem
.
integrate
(
gradu_vec
,
int_gradu_vec
,
spatial_dimension
*
spatial_dimension
,
element_type
,
_not_ghost
,
elem_filter
);
for
(
UInt
k
=
0
;
k
<
elem_filter
.
size
();
++
k
)
/// averaging is done only for normal components, so stress and strain
/// are equal
average
+=
0.5
*
(
int_gradu_vec
(
k
,
row_index
*
spatial_dimension
+
col_index
)
+
int_gradu_vec
(
k
,
col_index
*
spatial_dimension
+
row_index
));
}
}
else
if
(
field_type
==
"eigen_grad_u"
)
{
for
(
UInt
m
=
0
;
m
<
this
->
materials
.
size
();
++
m
)
{
const
auto
&
eigen_gradu_vec
=
this
->
materials
[
m
]
->
getInternal
<
Real
>
(
"eigen_grad_u"
)(
element_type
);
const
auto
&
elem_filter
=
this
->
materials
[
m
]
->
getElementFilter
(
element_type
);
Array
<
Real
>
int_eigen_gradu_vec
(
elem_filter
.
size
(),
spatial_dimension
*
spatial_dimension
,
"int_of_gradu"
);
fem
.
integrate
(
eigen_gradu_vec
,
int_eigen_gradu_vec
,
spatial_dimension
*
spatial_dimension
,
element_type
,
_not_ghost
,
elem_filter
);
for
(
UInt
k
=
0
;
k
<
elem_filter
.
size
();
++
k
)
/// averaging is done only for normal components, so stress and strain
/// are equal
average
+=
int_eigen_gradu_vec
(
k
,
row_index
*
spatial_dimension
+
col_index
);
}
}
else
{
AKANTU_ERROR
(
"Averaging not implemented for this field!!!"
);
}
}
return
average
/
this
->
volume
;
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModelRVE
::
homogenizeStiffness
(
Matrix
<
Real
>
&
C_macro
)
{
AKANTU_DEBUG_IN
();
const
Int
dim
=
2
;
AKANTU_DEBUG_ASSERT
(
this
->
spatial_dimension
==
dim
,
"Is only implemented for 2D!!!"
);
/// apply three independent loading states to determine C
/// 1. eps_el = (1;0;0) 2. eps_el = (0,1,0) 3. eps_el = (0,0,0.5)
/// clear the eigenstrain
Matrix
<
Real
>
zero_eigengradu
(
dim
,
dim
,
0.
);
GhostType
gt
=
_not_ghost
;
for
(
auto
element_type
:
mesh
.
elementTypes
(
dim
,
gt
,
_ek_not_defined
))
{
auto
&
prestrain_vect
=
const_cast
<
Array
<
Real
>
&>
(
this
->
getMaterial
(
"gel"
).
getInternal
<
Real
>
(
"eigen_grad_u"
)(
element_type
));
auto
prestrain_it
=
prestrain_vect
.
begin
(
spatial_dimension
,
spatial_dimension
);
auto
prestrain_end
=
prestrain_vect
.
end
(
spatial_dimension
,
spatial_dimension
);
for
(;
prestrain_it
!=
prestrain_end
;
++
prestrain_it
)
(
*
prestrain_it
)
=
zero_eigengradu
;
}
/// storage for results of 3 different loading states
UInt
voigt_size
=
VoigtHelper
<
dim
>::
size
;
Matrix
<
Real
>
stresses
(
voigt_size
,
voigt_size
,
0.
);
Matrix
<
Real
>
strains
(
voigt_size
,
voigt_size
,
0.
);
Matrix
<
Real
>
H
(
dim
,
dim
,
0.
);
/// save the damage state before filling up cracks
// ElementTypeMapReal saved_damage("saved_damage");
// saved_damage.initialize(getFEEngine(), _nb_component = 1, _default_value =
// 0);
// this->fillCracks(saved_damage);
/// virtual test 1:
H
(
0
,
0
)
=
0.01
;
this
->
performVirtualTesting
(
H
,
stresses
,
strains
,
0
);
/// virtual test 2:
H
.
zero
();
H
(
1
,
1
)
=
0.01
;
this
->
performVirtualTesting
(
H
,
stresses
,
strains
,
1
);
/// virtual test 3:
H
.
zero
();
H
(
0
,
1
)
=
0.01
;
this
->
performVirtualTesting
(
H
,
stresses
,
strains
,
2
);
/// drain cracks
// this->drainCracks(saved_damage);
/// compute effective stiffness
Matrix
<
Real
>
eps_inverse
(
voigt_size
,
voigt_size
);
eps_inverse
.
inverse
(
strains
);
C_macro
.
mul
<
false
,
false
>
(
stresses
,
eps_inverse
);
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModelRVE
::
performVirtualTesting
(
const
Matrix
<
Real
>
&
H
,
Matrix
<
Real
>
&
eff_stresses
,
Matrix
<
Real
>
&
eff_strains
,
const
UInt
test_no
)
{
AKANTU_DEBUG_IN
();
this
->
applyBoundaryConditions
(
H
);
auto
&
solver
=
this
->
getNonLinearSolver
();
solver
.
set
(
"max_iterations"
,
2
);
solver
.
set
(
"threshold"
,
1e-6
);
solver
.
set
(
"convergence_type"
,
SolveConvergenceCriteria
::
_solution
);
this
->
solveStep
();
/// get average stress and strain
eff_stresses
(
0
,
test_no
)
=
this
->
averageTensorField
(
0
,
0
,
"stress"
);
eff_strains
(
0
,
test_no
)
=
this
->
averageTensorField
(
0
,
0
,
"strain"
);
eff_stresses
(
1
,
test_no
)
=
this
->
averageTensorField
(
1
,
1
,
"stress"
);
eff_strains
(
1
,
test_no
)
=
this
->
averageTensorField
(
1
,
1
,
"strain"
);
eff_stresses
(
2
,
test_no
)
=
this
->
averageTensorField
(
1
,
0
,
"stress"
);
eff_strains
(
2
,
test_no
)
=
2.
*
this
->
averageTensorField
(
1
,
0
,
"strain"
);
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModelRVE
::
homogenizeEigenGradU
(
Matrix
<
Real
>
&
eigen_gradu_macro
)
{
AKANTU_DEBUG_IN
();
eigen_gradu_macro
(
0
,
0
)
=
this
->
averageTensorField
(
0
,
0
,
"eigen_grad_u"
);
eigen_gradu_macro
(
1
,
1
)
=
this
->
averageTensorField
(
1
,
1
,
"eigen_grad_u"
);
eigen_gradu_macro
(
0
,
1
)
=
this
->
averageTensorField
(
0
,
1
,
"eigen_grad_u"
);
eigen_gradu_macro
(
1
,
0
)
=
this
->
averageTensorField
(
1
,
0
,
"eigen_grad_u"
);
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModelRVE
::
initMaterials
()
{
AKANTU_DEBUG_IN
();
// make sure the material are instantiated
if
(
!
are_materials_instantiated
)
instantiateMaterials
();
if
(
use_RVE_mat_selector
)
{
const
Vector
<
Real
>
&
lowerBounds
=
mesh
.
getLowerBounds
();
const
Vector
<
Real
>
&
upperBounds
=
mesh
.
getUpperBounds
();
Real
bottom
=
lowerBounds
(
1
);
Real
top
=
upperBounds
(
1
);
Real
box_size
=
std
::
abs
(
top
-
bottom
);
Real
eps
=
box_size
*
1e-6
;
auto
tmp
=
std
::
make_shared
<
GelMaterialSelector
>
(
*
this
,
box_size
,
"gel"
,
this
->
nb_gel_pockets
,
eps
);
tmp
->
setFallback
(
material_selector
);
material_selector
=
tmp
;
}
this
->
assignMaterialToElements
();
// synchronize the element material arrays
this
->
synchronize
(
SynchronizationTag
::
_constitutive_law_id
);
for
(
auto
&
material
:
materials
)
{
/// init internals properties
const
auto
type
=
material
->
getID
();
if
(
type
.
find
(
"material_FE2"
)
!=
std
::
string
::
npos
)
continue
;
material
->
initMaterial
();
}
this
->
synchronize
(
SynchronizationTag
::
_clh_init_cl
);
if
(
this
->
non_local_manager
)
{
this
->
non_local_manager
->
initialize
();
}
// SolidMechanicsModel::initMaterials();
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModelRVE
::
fillCracks
(
ElementTypeMapReal
&
saved_damage
)
{
const
auto
&
mat_gel
=
this
->
getMaterial
(
"gel"
);
Real
E_gel
=
mat_gel
.
get
(
"E"
);
Real
E_homogenized
=
0.
;
for
(
auto
&&
mat
:
materials
)
{
if
(
mat
->
getName
()
==
"gel"
||
mat
->
getName
()
==
"FE2_mat"
)
continue
;
Real
E
=
mat
->
get
(
"E"
);
auto
&
damage
=
mat
->
getInternal
<
Real
>
(
"damage"
);
for
(
auto
&&
type
:
damage
.
elementTypes
())
{
const
auto
&
elem_filter
=
mat
->
getElementFilter
(
type
);
auto
nb_integration_point
=
getFEEngine
().
getNbIntegrationPoints
(
type
);
auto
sav_dam_it
=
make_view
(
saved_damage
(
type
),
nb_integration_point
).
begin
();
for
(
auto
&&
data
:
zip
(
elem_filter
,
make_view
(
damage
(
type
),
nb_integration_point
)))
{
auto
el
=
std
::
get
<
0
>
(
data
);
auto
&
dam
=
std
::
get
<
1
>
(
data
);
Vector
<
Real
>
sav_dam
=
sav_dam_it
[
el
];
sav_dam
=
dam
;
for
(
auto
q
:
arange
(
dam
.
size
()))
{
E_homogenized
=
(
E_gel
-
E
)
*
dam
(
q
)
+
E
;
dam
(
q
)
=
1.
-
(
E_homogenized
/
E
);
}
}
}
}
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModelRVE
::
drainCracks
(
const
ElementTypeMapReal
&
saved_damage
)
{
for
(
auto
&&
mat
:
materials
)
{
if
(
mat
->
getName
()
==
"gel"
||
mat
->
getName
()
==
"FE2_mat"
)
continue
;
auto
&
damage
=
mat
->
getInternal
<
Real
>
(
"damage"
);
for
(
auto
&&
type
:
damage
.
elementTypes
())
{
const
auto
&
elem_filter
=
mat
->
getElementFilter
(
type
);
auto
nb_integration_point
=
getFEEngine
().
getNbIntegrationPoints
(
type
);
auto
sav_dam_it
=
make_view
(
saved_damage
(
type
),
nb_integration_point
).
begin
();
for
(
auto
&&
data
:
zip
(
elem_filter
,
make_view
(
damage
(
type
),
nb_integration_point
)))
{
auto
el
=
std
::
get
<
0
>
(
data
);
auto
&
dam
=
std
::
get
<
1
>
(
data
);
Vector
<
Real
>
sav_dam
=
sav_dam_it
[
el
];
dam
=
sav_dam
;
}
}
}
}
}
// namespace akantu
Event Timeline
Log In to Comment