Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F79220471
material_reinforcement_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
Sun, Aug 25, 08:24
Size
26 KB
Mime Type
text/x-c++
Expires
Tue, Aug 27, 08:24 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
20161021
Attached To
rAKA akantu
material_reinforcement_tmpl.hh
View Options
/**
* @file material_reinforcement_tmpl.hh
*
* @author Lucas Frerot <lucas.frerot@epfl.ch>
*
* @date creation: Thu Feb 1 2018
*
* @brief Reinforcement material
*
* @section LICENSE
*
* Copyright (©) 2015 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 "aka_common.hh"
#include "aka_voigthelper.hh"
#include "material_reinforcement.hh"
namespace
akantu
{
/* -------------------------------------------------------------------------- */
template
<
class
Mat
,
UInt
dim
>
MaterialReinforcement
<
Mat
,
dim
>::
MaterialReinforcement
(
EmbeddedInterfaceModel
&
model
,
const
ID
&
id
)
:
Mat
(
model
,
1
,
model
.
getInterfaceMesh
(),
model
.
getFEEngine
(
"EmbeddedInterfaceFEEngine"
),
id
),
emodel
(
model
),
gradu_embedded
(
"gradu_embedded"
,
*
this
,
1
,
model
.
getFEEngine
(
"EmbeddedInterfaceFEEngine"
),
this
->
element_filter
),
directing_cosines
(
"directing_cosines"
,
*
this
,
1
,
model
.
getFEEngine
(
"EmbeddedInterfaceFEEngine"
),
this
->
element_filter
),
pre_stress
(
"pre_stress"
,
*
this
,
1
,
model
.
getFEEngine
(
"EmbeddedInterfaceFEEngine"
),
this
->
element_filter
),
area
(
1.0
),
shape_derivatives
()
{
AKANTU_DEBUG_IN
();
this
->
initialize
();
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
class
Mat
,
UInt
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
initialize
()
{
AKANTU_DEBUG_IN
();
this
->
registerParam
(
"area"
,
area
,
_pat_parsable
|
_pat_modifiable
,
"Reinforcement cross-sectional area"
);
this
->
registerParam
(
"pre_stress"
,
pre_stress
,
_pat_parsable
|
_pat_modifiable
,
"Uniform pre-stress"
);
// this->unregisterInternal(this->stress);
// Fool the AvgHomogenizingFunctor
// stress.initialize(dim * dim);
// Reallocate the element filter
this
->
element_filter
.
initialize
(
this
->
emodel
.
getInterfaceMesh
(),
_spatial_dimension
=
1
);
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
class
Mat
,
UInt
dim
>
MaterialReinforcement
<
Mat
,
dim
>::~
MaterialReinforcement
()
{
AKANTU_DEBUG_IN
();
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
class
Mat
,
UInt
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
initMaterial
()
{
Mat
::
initMaterial
();
gradu_embedded
.
initialize
(
dim
*
dim
);
pre_stress
.
initialize
(
1
);
/// We initialise the stuff that is not going to change during the simulation
this
->
initFilters
();
this
->
allocBackgroundShapeDerivatives
();
this
->
initBackgroundShapeDerivatives
();
this
->
initDirectingCosines
();
}
/* -------------------------------------------------------------------------- */
/// Initialize the filter for background elements
template
<
class
Mat
,
UInt
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
initFilters
()
{
for
(
auto
gt
:
ghost_types
)
{
for
(
auto
&&
type
:
emodel
.
getInterfaceMesh
().
elementTypes
(
1
,
gt
))
{
std
::
string
shaped_id
=
"filter"
;
if
(
gt
==
_ghost
)
shaped_id
+=
":ghost"
;
auto
&
background
=
background_filter
(
std
::
make_unique
<
ElementTypeMapArray
<
UInt
>>
(
"bg_"
+
shaped_id
,
this
->
name
),
type
,
gt
);
auto
&
foreground
=
foreground_filter
(
std
::
make_unique
<
ElementTypeMapArray
<
UInt
>>
(
shaped_id
,
this
->
name
),
type
,
gt
);
foreground
->
initialize
(
emodel
.
getMesh
(),
_nb_component
=
1
,
_ghost_type
=
gt
);
background
->
initialize
(
emodel
.
getMesh
(),
_nb_component
=
1
,
_ghost_type
=
gt
);
// Computing filters
for
(
auto
&&
bg_type
:
background
->
elementTypes
(
dim
,
gt
))
{
filterInterfaceBackgroundElements
(
(
*
foreground
)(
bg_type
),
(
*
background
)(
bg_type
),
bg_type
,
type
,
gt
);
}
}
}
}
/* -------------------------------------------------------------------------- */
/// Construct a filter for a (interface_type, background_type) pair
template
<
class
Mat
,
UInt
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
filterInterfaceBackgroundElements
(
Array
<
UInt
>
&
foreground
,
Array
<
UInt
>
&
background
,
const
ElementType
&
type
,
const
ElementType
&
interface_type
,
GhostType
ghost_type
)
{
AKANTU_DEBUG_IN
();
foreground
.
resize
(
0
);
background
.
resize
(
0
);
Array
<
Element
>
&
elements
=
emodel
.
getInterfaceAssociatedElements
(
interface_type
,
ghost_type
);
Array
<
UInt
>
&
elem_filter
=
this
->
element_filter
(
interface_type
,
ghost_type
);
for
(
auto
&
elem_id
:
elem_filter
)
{
Element
&
elem
=
elements
(
elem_id
);
if
(
elem
.
type
==
type
)
{
background
.
push_back
(
elem
.
element
);
foreground
.
push_back
(
elem_id
);
}
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
namespace
detail
{
class
BackgroundShapeDInitializer
:
public
ElementTypeMapArrayInitializer
{
public
:
BackgroundShapeDInitializer
(
UInt
spatial_dimension
,
FEEngine
&
engine
,
const
ElementType
&
foreground_type
,
const
ElementTypeMapArray
<
UInt
>
&
filter
,
const
GhostType
&
ghost_type
)
:
ElementTypeMapArrayInitializer
(
[](
const
ElementType
&
bgtype
,
const
GhostType
&
)
{
return
ShapeFunctions
::
getShapeDerivativesSize
(
bgtype
);
},
spatial_dimension
,
ghost_type
,
_ek_regular
)
{
auto
nb_quad
=
engine
.
getNbIntegrationPoints
(
foreground_type
);
// Counting how many background elements are affected by elements of
// interface_type
for
(
auto
type
:
filter
.
elementTypes
(
this
->
spatial_dimension
))
{
// Inserting size
array_size_per_bg_type
(
filter
(
type
).
size
()
*
nb_quad
,
type
,
this
->
ghost_type
);
}
}
auto
elementTypes
()
const
->
decltype
(
auto
)
{
return
array_size_per_bg_type
.
elementTypes
();
}
UInt
size
(
const
ElementType
&
bgtype
)
const
{
return
array_size_per_bg_type
(
bgtype
,
this
->
ghost_type
);
}
protected
:
ElementTypeMap
<
UInt
>
array_size_per_bg_type
;
};
}
/**
* Background shape derivatives need to be stored per background element
* types but also per embedded element type, which is why they are stored
* in an ElementTypeMap<ElementTypeMapArray<Real> *>. The outer ElementTypeMap
* refers to the embedded types, and the inner refers to the background types.
*/
template
<
class
Mat
,
UInt
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
allocBackgroundShapeDerivatives
()
{
AKANTU_DEBUG_IN
();
for
(
auto
gt
:
ghost_types
)
{
for
(
auto
&&
type
:
emodel
.
getInterfaceMesh
().
elementTypes
(
1
,
gt
))
{
std
::
string
shaped_id
=
"embedded_shape_derivatives"
;
if
(
gt
==
_ghost
)
shaped_id
+=
":ghost"
;
auto
&
shaped_etma
=
shape_derivatives
(
std
::
make_unique
<
ElementTypeMapArray
<
Real
>>
(
shaped_id
,
this
->
name
),
type
,
gt
);
shaped_etma
->
initialize
(
detail
::
BackgroundShapeDInitializer
(
emodel
.
getSpatialDimension
(),
emodel
.
getFEEngine
(
"EmbeddedInterfaceFEEngine"
),
type
,
*
background_filter
(
type
,
gt
),
gt
),
0
,
true
);
}
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
class
Mat
,
UInt
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
initBackgroundShapeDerivatives
()
{
AKANTU_DEBUG_IN
();
for
(
auto
interface_type
:
this
->
element_filter
.
elementTypes
(
this
->
spatial_dimension
))
{
for
(
auto
type
:
background_filter
(
interface_type
)
->
elementTypes
(
dim
))
{
computeBackgroundShapeDerivatives
(
interface_type
,
type
,
_not_ghost
,
this
->
element_filter
(
interface_type
));
}
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
class
Mat
,
UInt
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
computeBackgroundShapeDerivatives
(
const
ElementType
&
interface_type
,
const
ElementType
&
bg_type
,
GhostType
ghost_type
,
const
Array
<
UInt
>
&
filter
)
{
auto
&
interface_engine
=
emodel
.
getFEEngine
(
"EmbeddedInterfaceFEEngine"
);
auto
&
engine
=
emodel
.
getFEEngine
();
auto
&
interface_mesh
=
emodel
.
getInterfaceMesh
();
const
auto
nb_nodes_elem_bg
=
Mesh
::
getNbNodesPerElement
(
bg_type
);
// const auto nb_strss = VoigtHelper<dim>::size;
const
auto
nb_quads_per_elem
=
interface_engine
.
getNbIntegrationPoints
(
interface_type
);
Array
<
Real
>
quad_pos
(
0
,
dim
,
"interface_quad_pos"
);
interface_engine
.
interpolateOnIntegrationPoints
(
interface_mesh
.
getNodes
(),
quad_pos
,
dim
,
interface_type
,
ghost_type
,
filter
);
auto
&
background_shapesd
=
(
*
shape_derivatives
(
interface_type
,
ghost_type
))(
bg_type
,
ghost_type
);
auto
&
background_elements
=
(
*
background_filter
(
interface_type
,
ghost_type
))(
bg_type
,
ghost_type
);
auto
&
foreground_elements
=
(
*
foreground_filter
(
interface_type
,
ghost_type
))(
bg_type
,
ghost_type
);
auto
shapesd_begin
=
background_shapesd
.
begin
(
dim
,
nb_nodes_elem_bg
,
nb_quads_per_elem
);
auto
quad_begin
=
quad_pos
.
begin
(
dim
,
nb_quads_per_elem
);
for
(
auto
&&
tuple
:
zip
(
background_elements
,
foreground_elements
))
{
UInt
bg
=
std
::
get
<
0
>
(
tuple
),
fg
=
std
::
get
<
1
>
(
tuple
);
for
(
UInt
i
=
0
;
i
<
nb_quads_per_elem
;
++
i
)
{
Matrix
<
Real
>
shapesd
=
Tensor3
<
Real
>
(
shapesd_begin
[
fg
])(
i
);
Vector
<
Real
>
quads
=
Matrix
<
Real
>
(
quad_begin
[
fg
])(
i
);
engine
.
computeShapeDerivatives
(
quads
,
bg
,
bg_type
,
shapesd
,
ghost_type
);
}
}
}
/* -------------------------------------------------------------------------- */
template
<
class
Mat
,
UInt
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
initDirectingCosines
()
{
AKANTU_DEBUG_IN
();
Mesh
&
mesh
=
emodel
.
getInterfaceMesh
();
Mesh
::
type_iterator
type_it
=
mesh
.
firstType
(
1
,
_not_ghost
);
Mesh
::
type_iterator
type_end
=
mesh
.
lastType
(
1
,
_not_ghost
);
const
UInt
voigt_size
=
VoigtHelper
<
dim
>::
size
;
directing_cosines
.
initialize
(
voigt_size
);
for
(;
type_it
!=
type_end
;
++
type_it
)
{
computeDirectingCosines
(
*
type_it
,
_not_ghost
);
// computeDirectingCosines(*type_it, _ghost);
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
class
Mat
,
UInt
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
assembleStiffnessMatrix
(
GhostType
ghost_type
)
{
AKANTU_DEBUG_IN
();
Mesh
&
interface_mesh
=
emodel
.
getInterfaceMesh
();
Mesh
::
type_iterator
type_it
=
interface_mesh
.
firstType
(
1
,
_not_ghost
);
Mesh
::
type_iterator
type_end
=
interface_mesh
.
lastType
(
1
,
_not_ghost
);
for
(;
type_it
!=
type_end
;
++
type_it
)
{
assembleStiffnessMatrix
(
*
type_it
,
ghost_type
);
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
class
Mat
,
UInt
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
assembleInternalForces
(
GhostType
ghost_type
)
{
AKANTU_DEBUG_IN
();
Mesh
&
interface_mesh
=
emodel
.
getInterfaceMesh
();
Mesh
::
type_iterator
type_it
=
interface_mesh
.
firstType
(
1
,
_not_ghost
);
Mesh
::
type_iterator
type_end
=
interface_mesh
.
lastType
(
1
,
_not_ghost
);
for
(;
type_it
!=
type_end
;
++
type_it
)
{
this
->
assembleInternalForces
(
*
type_it
,
ghost_type
);
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
class
Mat
,
UInt
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
computeAllStresses
(
GhostType
ghost_type
)
{
AKANTU_DEBUG_IN
();
Mesh
::
type_iterator
it
=
emodel
.
getInterfaceMesh
().
firstType
();
Mesh
::
type_iterator
last_type
=
emodel
.
getInterfaceMesh
().
lastType
();
for
(;
it
!=
last_type
;
++
it
)
{
computeGradU
(
*
it
,
ghost_type
);
this
->
computeStress
(
*
it
,
ghost_type
);
addPrestress
(
*
it
,
ghost_type
);
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
class
Mat
,
UInt
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
addPrestress
(
const
ElementType
&
type
,
GhostType
ghost_type
)
{
auto
&
stress
=
this
->
stress
(
type
,
ghost_type
);
auto
&
sigma_p
=
this
->
pre_stress
(
type
,
ghost_type
);
for
(
auto
&&
tuple
:
zip
(
stress
,
sigma_p
))
{
std
::
get
<
0
>
(
tuple
)
+=
std
::
get
<
1
>
(
tuple
);
}
}
/* -------------------------------------------------------------------------- */
template
<
class
Mat
,
UInt
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
assembleInternalForces
(
const
ElementType
&
type
,
GhostType
ghost_type
)
{
AKANTU_DEBUG_IN
();
Mesh
&
mesh
=
emodel
.
getMesh
();
Mesh
::
type_iterator
type_it
=
mesh
.
firstType
(
dim
,
ghost_type
);
Mesh
::
type_iterator
type_end
=
mesh
.
lastType
(
dim
,
ghost_type
);
for
(;
type_it
!=
type_end
;
++
type_it
)
{
assembleInternalForcesInterface
(
type
,
*
type_it
,
ghost_type
);
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
/**
* Computes and assemble the residual. Residual in reinforcement is computed as:
*
* \f[
* \vec{r} = A_s \int_S{\mathbf{B}^T\mathbf{C}^T \vec{\sigma_s}\,\mathrm{d}s}
* \f]
*/
template
<
class
Mat
,
UInt
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
assembleInternalForcesInterface
(
const
ElementType
&
interface_type
,
const
ElementType
&
background_type
,
GhostType
ghost_type
)
{
AKANTU_DEBUG_IN
();
UInt
voigt_size
=
VoigtHelper
<
dim
>::
size
;
FEEngine
&
interface_engine
=
emodel
.
getFEEngine
(
"EmbeddedInterfaceFEEngine"
);
Array
<
UInt
>
&
elem_filter
=
this
->
element_filter
(
interface_type
,
ghost_type
);
UInt
nodes_per_background_e
=
Mesh
::
getNbNodesPerElement
(
background_type
);
UInt
nb_quadrature_points
=
interface_engine
.
getNbIntegrationPoints
(
interface_type
,
ghost_type
);
UInt
nb_element
=
elem_filter
.
size
();
UInt
back_dof
=
dim
*
nodes_per_background_e
;
Array
<
Real
>
&
shapesd
=
(
*
shape_derivatives
(
interface_type
,
ghost_type
))(
background_type
,
ghost_type
);
Array
<
Real
>
integrant
(
nb_quadrature_points
*
nb_element
,
back_dof
,
"integrant"
);
auto
integrant_it
=
integrant
.
begin
(
back_dof
);
auto
integrant_end
=
integrant
.
end
(
back_dof
);
Array
<
Real
>::
matrix_iterator
B_it
=
shapesd
.
begin
(
dim
,
nodes_per_background_e
);
auto
C_it
=
directing_cosines
(
interface_type
,
ghost_type
).
begin
(
voigt_size
);
auto
sigma_it
=
this
->
stress
(
interface_type
,
ghost_type
).
begin
();
Matrix
<
Real
>
Bvoigt
(
voigt_size
,
back_dof
);
for
(;
integrant_it
!=
integrant_end
;
++
integrant_it
,
++
B_it
,
++
C_it
,
++
sigma_it
)
{
VoigtHelper
<
dim
>::
transferBMatrixToSymVoigtBMatrix
(
*
B_it
,
Bvoigt
,
nodes_per_background_e
);
Vector
<
Real
>
&
C
=
*
C_it
;
Vector
<
Real
>
&
BtCt_sigma
=
*
integrant_it
;
BtCt_sigma
.
mul
<
true
>
(
Bvoigt
,
C
);
BtCt_sigma
*=
*
sigma_it
*
area
;
}
Array
<
Real
>
residual_interface
(
nb_element
,
back_dof
,
"residual_interface"
);
interface_engine
.
integrate
(
integrant
,
residual_interface
,
back_dof
,
interface_type
,
ghost_type
,
elem_filter
);
integrant
.
resize
(
0
);
Array
<
UInt
>
background_filter
(
nb_element
,
1
,
"background_filter"
);
auto
&
filter
=
getBackgroundFilter
(
interface_type
,
background_type
,
ghost_type
);
emodel
.
getDOFManager
().
assembleElementalArrayLocalArray
(
residual_interface
,
emodel
.
getInternalForce
(),
background_type
,
ghost_type
,
-
1.
,
filter
);
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
class
Mat
,
UInt
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
computeDirectingCosines
(
const
ElementType
&
type
,
GhostType
ghost_type
)
{
AKANTU_DEBUG_IN
();
Mesh
&
interface_mesh
=
emodel
.
getInterfaceMesh
();
const
UInt
nb_nodes_per_element
=
Mesh
::
getNbNodesPerElement
(
type
);
const
UInt
steel_dof
=
dim
*
nb_nodes_per_element
;
const
UInt
voigt_size
=
VoigtHelper
<
dim
>::
size
;
const
UInt
nb_quad_points
=
emodel
.
getFEEngine
(
"EmbeddedInterfaceFEEngine"
)
.
getNbIntegrationPoints
(
type
,
ghost_type
);
Array
<
Real
>
node_coordinates
(
this
->
element_filter
(
type
,
ghost_type
).
size
(),
steel_dof
);
this
->
emodel
.
getFEEngine
().
template
extractNodalToElementField
<
Real
>
(
interface_mesh
,
interface_mesh
.
getNodes
(),
node_coordinates
,
type
,
ghost_type
,
this
->
element_filter
(
type
,
ghost_type
));
Array
<
Real
>::
matrix_iterator
directing_cosines_it
=
directing_cosines
(
type
,
ghost_type
).
begin
(
1
,
voigt_size
);
Array
<
Real
>::
matrix_iterator
node_coordinates_it
=
node_coordinates
.
begin
(
dim
,
nb_nodes_per_element
);
Array
<
Real
>::
matrix_iterator
node_coordinates_end
=
node_coordinates
.
end
(
dim
,
nb_nodes_per_element
);
for
(;
node_coordinates_it
!=
node_coordinates_end
;
++
node_coordinates_it
)
{
for
(
UInt
i
=
0
;
i
<
nb_quad_points
;
i
++
,
++
directing_cosines_it
)
{
Matrix
<
Real
>
&
nodes
=
*
node_coordinates_it
;
Matrix
<
Real
>
&
cosines
=
*
directing_cosines_it
;
computeDirectingCosinesOnQuad
(
nodes
,
cosines
);
}
}
// Mauro: the directing_cosines internal is defined on the quadrature points
// of each element
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
class
Mat
,
UInt
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
assembleStiffnessMatrix
(
const
ElementType
&
type
,
GhostType
ghost_type
)
{
AKANTU_DEBUG_IN
();
Mesh
&
mesh
=
emodel
.
getMesh
();
Mesh
::
type_iterator
type_it
=
mesh
.
firstType
(
dim
,
ghost_type
);
Mesh
::
type_iterator
type_end
=
mesh
.
lastType
(
dim
,
ghost_type
);
for
(;
type_it
!=
type_end
;
++
type_it
)
{
assembleStiffnessMatrixInterface
(
type
,
*
type_it
,
ghost_type
);
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
/**
* Computes the reinforcement stiffness matrix (Gomes & Awruch, 2001)
* \f[
* \mathbf{K}_e = \sum_{i=1}^R{A_i\int_{S_i}{\mathbf{B}^T
* \mathbf{C}_i^T \mathbf{D}_{s, i} \mathbf{C}_i \mathbf{B}\,\mathrm{d}s}}
* \f]
*/
template
<
class
Mat
,
UInt
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
assembleStiffnessMatrixInterface
(
const
ElementType
&
interface_type
,
const
ElementType
&
background_type
,
GhostType
ghost_type
)
{
AKANTU_DEBUG_IN
();
UInt
voigt_size
=
VoigtHelper
<
dim
>::
size
;
FEEngine
&
interface_engine
=
emodel
.
getFEEngine
(
"EmbeddedInterfaceFEEngine"
);
Array
<
UInt
>
&
elem_filter
=
this
->
element_filter
(
interface_type
,
ghost_type
);
Array
<
Real
>
&
grad_u
=
gradu_embedded
(
interface_type
,
ghost_type
);
UInt
nb_element
=
elem_filter
.
size
();
UInt
nodes_per_background_e
=
Mesh
::
getNbNodesPerElement
(
background_type
);
UInt
nb_quadrature_points
=
interface_engine
.
getNbIntegrationPoints
(
interface_type
,
ghost_type
);
UInt
back_dof
=
dim
*
nodes_per_background_e
;
UInt
integrant_size
=
back_dof
;
grad_u
.
resize
(
nb_quadrature_points
*
nb_element
);
Array
<
Real
>
tangent_moduli
(
nb_element
*
nb_quadrature_points
,
1
,
"interface_tangent_moduli"
);
this
->
computeTangentModuli
(
interface_type
,
tangent_moduli
,
ghost_type
);
Array
<
Real
>
&
shapesd
=
(
*
shape_derivatives
(
interface_type
,
ghost_type
))(
background_type
,
ghost_type
);
Array
<
Real
>
integrant
(
nb_element
*
nb_quadrature_points
,
integrant_size
*
integrant_size
,
"B^t*C^t*D*C*B"
);
/// Temporary matrices for integrant product
Matrix
<
Real
>
Bvoigt
(
voigt_size
,
back_dof
);
Matrix
<
Real
>
DCB
(
1
,
back_dof
);
Matrix
<
Real
>
CtDCB
(
voigt_size
,
back_dof
);
Array
<
Real
>::
scalar_iterator
D_it
=
tangent_moduli
.
begin
();
Array
<
Real
>::
scalar_iterator
D_end
=
tangent_moduli
.
end
();
Array
<
Real
>::
matrix_iterator
C_it
=
directing_cosines
(
interface_type
,
ghost_type
).
begin
(
1
,
voigt_size
);
Array
<
Real
>::
matrix_iterator
B_it
=
shapesd
.
begin
(
dim
,
nodes_per_background_e
);
Array
<
Real
>::
matrix_iterator
integrant_it
=
integrant
.
begin
(
integrant_size
,
integrant_size
);
for
(;
D_it
!=
D_end
;
++
D_it
,
++
C_it
,
++
B_it
,
++
integrant_it
)
{
Real
&
D
=
*
D_it
;
Matrix
<
Real
>
&
C
=
*
C_it
;
Matrix
<
Real
>
&
B
=
*
B_it
;
Matrix
<
Real
>
&
BtCtDCB
=
*
integrant_it
;
VoigtHelper
<
dim
>::
transferBMatrixToSymVoigtBMatrix
(
B
,
Bvoigt
,
nodes_per_background_e
);
DCB
.
mul
<
false
,
false
>
(
C
,
Bvoigt
);
DCB
*=
D
*
area
;
CtDCB
.
mul
<
true
,
false
>
(
C
,
DCB
);
BtCtDCB
.
mul
<
true
,
false
>
(
Bvoigt
,
CtDCB
);
}
tangent_moduli
.
resize
(
0
);
Array
<
Real
>
K_interface
(
nb_element
,
integrant_size
*
integrant_size
,
"K_interface"
);
interface_engine
.
integrate
(
integrant
,
K_interface
,
integrant_size
*
integrant_size
,
interface_type
,
ghost_type
,
elem_filter
);
integrant
.
resize
(
0
);
// Mauro: Here K_interface contains the local stiffness matrices,
// directing_cosines contains the information about the orientation
// of the reinforcements, any rotation of the local stiffness matrix
// can be done here
auto
&
filter
=
getBackgroundFilter
(
interface_type
,
background_type
,
ghost_type
);
emodel
.
getDOFManager
().
assembleElementalMatricesToMatrix
(
"K"
,
"displacement"
,
K_interface
,
background_type
,
ghost_type
,
_symmetric
,
filter
);
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
class
Mat
,
UInt
dim
>
Real
MaterialReinforcement
<
Mat
,
dim
>::
getEnergy
(
const
std
::
string
&
id
)
{
AKANTU_DEBUG_IN
();
if
(
id
==
"potential"
)
{
Real
epot
=
0.
;
this
->
computePotentialEnergyByElements
();
Mesh
::
type_iterator
it
=
this
->
element_filter
.
firstType
(
this
->
spatial_dimension
),
end
=
this
->
element_filter
.
lastType
(
this
->
spatial_dimension
);
for
(;
it
!=
end
;
++
it
)
{
FEEngine
&
interface_engine
=
emodel
.
getFEEngine
(
"EmbeddedInterfaceFEEngine"
);
epot
+=
interface_engine
.
integrate
(
this
->
potential_energy
(
*
it
,
_not_ghost
),
*
it
,
_not_ghost
,
this
->
element_filter
(
*
it
,
_not_ghost
));
epot
*=
area
;
}
return
epot
;
}
AKANTU_DEBUG_OUT
();
return
0
;
}
/* -------------------------------------------------------------------------- */
template
<
class
Mat
,
UInt
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
computeGradU
(
const
ElementType
&
interface_type
,
GhostType
ghost_type
)
{
// Looping over background types
for
(
auto
&&
bg_type
:
background_filter
(
interface_type
,
ghost_type
)
->
elementTypes
(
dim
))
{
const
UInt
nodes_per_background_e
=
Mesh
::
getNbNodesPerElement
(
bg_type
);
const
UInt
voigt_size
=
VoigtHelper
<
dim
>::
size
;
auto
&
bg_shapesd
=
(
*
shape_derivatives
(
interface_type
,
ghost_type
))(
bg_type
,
ghost_type
);
auto
&
filter
=
getBackgroundFilter
(
interface_type
,
bg_type
,
ghost_type
);
Array
<
Real
>
disp_per_element
(
0
,
dim
*
nodes_per_background_e
,
"disp_elem"
);
FEEngine
::
extractNodalToElementField
(
emodel
.
getMesh
(),
emodel
.
getDisplacement
(),
disp_per_element
,
bg_type
,
ghost_type
,
filter
);
Matrix
<
Real
>
concrete_du
(
dim
,
dim
);
Matrix
<
Real
>
epsilon
(
dim
,
dim
);
Vector
<
Real
>
evoigt
(
voigt_size
);
for
(
auto
&&
tuple
:
zip
(
make_view
(
disp_per_element
,
dim
,
nodes_per_background_e
),
make_view
(
bg_shapesd
,
dim
,
nodes_per_background_e
),
this
->
gradu
(
interface_type
,
ghost_type
),
make_view
(
this
->
directing_cosines
(
interface_type
,
ghost_type
),
voigt_size
)))
{
auto
&
u
=
std
::
get
<
0
>
(
tuple
);
auto
&
B
=
std
::
get
<
1
>
(
tuple
);
auto
&
du
=
std
::
get
<
2
>
(
tuple
);
auto
&
C
=
std
::
get
<
3
>
(
tuple
);
concrete_du
.
mul
<
false
,
true
>
(
u
,
B
);
auto
epsilon
=
0.5
*
(
concrete_du
+
concrete_du
.
transpose
());
strainTensorToVoigtVector
(
epsilon
,
evoigt
);
du
=
C
.
dot
(
evoigt
);
}
}
}
/* -------------------------------------------------------------------------- */
/**
* The structure of the directing cosines matrix is :
* \f{eqnarray*}{
* C_{1,\cdot} & = & (l^2, m^2, n^2, mn, ln, lm) \\
* C_{i,j} & = & 0
* \f}
*
* with :
* \f[
* (l, m, n) = \frac{1}{\|\frac{\mathrm{d}\vec{r}(s)}{\mathrm{d}s}\|} \cdot
* \frac{\mathrm{d}\vec{r}(s)}{\mathrm{d}s}
* \f]
*/
template
<
class
Mat
,
UInt
dim
>
inline
void
MaterialReinforcement
<
Mat
,
dim
>::
computeDirectingCosinesOnQuad
(
const
Matrix
<
Real
>
&
nodes
,
Matrix
<
Real
>
&
cosines
)
{
AKANTU_DEBUG_IN
();
AKANTU_DEBUG_ASSERT
(
nodes
.
cols
()
==
2
,
"Higher order reinforcement elements not implemented"
);
const
Vector
<
Real
>
a
=
nodes
(
0
),
b
=
nodes
(
1
);
Vector
<
Real
>
delta
=
b
-
a
;
Real
sq_length
=
0.
;
for
(
UInt
i
=
0
;
i
<
dim
;
i
++
)
{
sq_length
+=
delta
(
i
)
*
delta
(
i
);
}
if
(
dim
==
2
)
{
cosines
(
0
,
0
)
=
delta
(
0
)
*
delta
(
0
);
// l^2
cosines
(
0
,
1
)
=
delta
(
1
)
*
delta
(
1
);
// m^2
cosines
(
0
,
2
)
=
delta
(
0
)
*
delta
(
1
);
// lm
}
else
if
(
dim
==
3
)
{
cosines
(
0
,
0
)
=
delta
(
0
)
*
delta
(
0
);
// l^2
cosines
(
0
,
1
)
=
delta
(
1
)
*
delta
(
1
);
// m^2
cosines
(
0
,
2
)
=
delta
(
2
)
*
delta
(
2
);
// n^2
cosines
(
0
,
3
)
=
delta
(
1
)
*
delta
(
2
);
// mn
cosines
(
0
,
4
)
=
delta
(
0
)
*
delta
(
2
);
// ln
cosines
(
0
,
5
)
=
delta
(
0
)
*
delta
(
1
);
// lm
}
cosines
/=
sq_length
;
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
class
Mat
,
UInt
dim
>
inline
void
MaterialReinforcement
<
Mat
,
dim
>::
stressTensorToVoigtVector
(
const
Matrix
<
Real
>
&
tensor
,
Vector
<
Real
>
&
vector
)
{
AKANTU_DEBUG_IN
();
for
(
UInt
i
=
0
;
i
<
dim
;
i
++
)
{
vector
(
i
)
=
tensor
(
i
,
i
);
}
if
(
dim
==
2
)
{
vector
(
2
)
=
tensor
(
0
,
1
);
}
else
if
(
dim
==
3
)
{
vector
(
3
)
=
tensor
(
1
,
2
);
vector
(
4
)
=
tensor
(
0
,
2
);
vector
(
5
)
=
tensor
(
0
,
1
);
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
class
Mat
,
UInt
dim
>
inline
void
MaterialReinforcement
<
Mat
,
dim
>::
strainTensorToVoigtVector
(
const
Matrix
<
Real
>
&
tensor
,
Vector
<
Real
>
&
vector
)
{
AKANTU_DEBUG_IN
();
for
(
UInt
i
=
0
;
i
<
dim
;
i
++
)
{
vector
(
i
)
=
tensor
(
i
,
i
);
}
if
(
dim
==
2
)
{
vector
(
2
)
=
2
*
tensor
(
0
,
1
);
}
else
if
(
dim
==
3
)
{
vector
(
3
)
=
2
*
tensor
(
1
,
2
);
vector
(
4
)
=
2
*
tensor
(
0
,
2
);
vector
(
5
)
=
2
*
tensor
(
0
,
1
);
}
AKANTU_DEBUG_OUT
();
}
}
// akantu
Event Timeline
Log In to Comment