Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F82446634
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
Wed, Sep 11, 11:29
Size
22 KB
Mime Type
text/x-c++
Expires
Fri, Sep 13, 11:29 (2 d)
Engine
blob
Format
Raw Data
Handle
20702021
Attached To
rAKA akantu
material_reinforcement_tmpl.hh
View Options
/**
* Copyright (©) 2015-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 "aka_common.hh"
#include "aka_voigthelper.hh"
#include "material_reinforcement.hh"
namespace
akantu
{
/* -------------------------------------------------------------------------- */
template
<
class
Mat
,
Int
dim
>
MaterialReinforcement
<
Mat
,
dim
>::
MaterialReinforcement
(
EmbeddedInterfaceModel
&
model
,
const
ID
&
id
)
:
Mat
(
model
,
id
,
"EmbeddedInterfaceFEEngine"
),
emodel
(
model
),
area
(
1.0
)
{
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"
);
// Reallocate the element filter
this
->
element_filter
.
initialize
(
this
->
emodel
.
getInterfaceMesh
(),
_spatial_dimension
=
1
);
}
/* -------------------------------------------------------------------------- */
template
<
class
Mat
,
Int
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
initMaterial
()
{
Mat
::
initMaterial
();
gradu_embedded
=
this
->
registerInternal
(
"gradu_embedded"
,
dim
*
dim
,
"EmbeddedInterfaceFEEngine"
);
pre_stress
=
this
->
registerInternal
(
"pre_stress"
,
1
,
"EmbeddedInterfaceFEEngine"
);
/// 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
,
Int
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
<
Idx
>>
(
"bg_"
+
shaped_id
,
this
->
name
),
type
,
gt
);
auto
&
foreground
=
foreground_filter
(
std
::
make_unique
<
ElementTypeMapArray
<
Idx
>>
(
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
,
Int
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
filterInterfaceBackgroundElements
(
Array
<
Idx
>
&
foreground
,
Array
<
Idx
>
&
background
,
ElementType
type
,
ElementType
interface_type
,
GhostType
ghost_type
)
{
AKANTU_DEBUG_IN
();
foreground
.
resize
(
0
);
background
.
resize
(
0
);
Array
<
Element
>
&
elements
=
emodel
.
getInterfaceAssociatedElements
(
interface_type
,
ghost_type
);
auto
&
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
(
Int
spatial_dimension
,
FEEngine
&
engine
,
ElementType
foreground_type
,
const
ElementTypeMapArray
<
Idx
>
&
filter
,
GhostType
ghost_type
)
:
ElementTypeMapArrayInitializer
(
[](
ElementType
bgtype
,
GhostType
/*unused*/
)
{
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
(
ElementType
bgtype
)
const
{
return
array_size_per_bg_type
(
bgtype
,
this
->
ghost_type
);
}
protected
:
ElementTypeMap
<
UInt
>
array_size_per_bg_type
;
};
}
// namespace detail
/**
* 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
,
Int
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
,
Int
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
,
Int
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
computeBackgroundShapeDerivatives
(
ElementType
interface_type
,
ElementType
bg_type
,
GhostType
ghost_type
,
const
Array
<
Idx
>
&
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
))
{
auto
bg
=
std
::
get
<
0
>
(
tuple
);
auto
fg
=
std
::
get
<
1
>
(
tuple
);
for
(
Int
i
=
0
;
i
<
nb_quads_per_elem
;
++
i
)
{
engine
.
computeShapeDerivatives
(
quad_begin
[
fg
](
i
),
bg
,
bg_type
,
shapesd_begin
[
fg
](
i
),
ghost_type
);
}
}
}
/* -------------------------------------------------------------------------- */
template
<
class
Mat
,
Int
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
initDirectingCosines
()
{
AKANTU_DEBUG_IN
();
Mesh
&
mesh
=
emodel
.
getInterfaceMesh
();
const
auto
voigt_size
=
VoigtHelper
<
dim
>::
size
;
directing_cosines
=
this
->
registerInternal
(
"directing_cosines"
,
voigt_size
,
"EmbeddedInterfaceFEEngine"
);
for
(
auto
&&
type
:
mesh
.
elementTypes
(
1
,
_not_ghost
))
{
computeDirectingCosines
(
type
,
_not_ghost
);
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
class
Mat
,
Int
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
assembleStiffnessMatrix
(
GhostType
ghost_type
)
{
AKANTU_DEBUG_IN
();
Mesh
&
interface_mesh
=
emodel
.
getInterfaceMesh
();
for
(
auto
&&
type
:
interface_mesh
.
elementTypes
(
1
,
_not_ghost
))
{
assembleStiffnessMatrix
(
type
,
ghost_type
);
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
class
Mat
,
Int
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
assembleInternalForces
(
GhostType
ghost_type
)
{
AKANTU_DEBUG_IN
();
Mesh
&
interface_mesh
=
emodel
.
getInterfaceMesh
();
for
(
auto
&&
type
:
interface_mesh
.
elementTypes
(
1
,
_not_ghost
))
{
this
->
assembleInternalForces
(
type
,
ghost_type
);
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
class
Mat
,
Int
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
computeAllStresses
(
GhostType
ghost_type
)
{
AKANTU_DEBUG_IN
();
Mesh
&
interface_mesh
=
emodel
.
getInterfaceMesh
();
for
(
auto
&&
type
:
interface_mesh
.
elementTypes
(
_ghost_type
=
ghost_type
))
{
computeGradU
(
type
,
ghost_type
);
this
->
computeStress
(
type
,
ghost_type
);
addPrestress
(
type
,
ghost_type
);
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
class
Mat
,
Int
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
addPrestress
(
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
,
Int
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
assembleInternalForces
(
ElementType
type
,
GhostType
ghost_type
)
{
AKANTU_DEBUG_IN
();
Mesh
&
mesh
=
emodel
.
getMesh
();
for
(
auto
&&
mesh_type
:
mesh
.
elementTypes
(
dim
,
ghost_type
))
{
assembleInternalForcesInterface
(
type
,
mesh_type
,
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
,
Int
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
assembleInternalForcesInterface
(
ElementType
interface_type
,
ElementType
background_type
,
GhostType
ghost_type
)
{
AKANTU_DEBUG_IN
();
constexpr
auto
voigt_size
=
VoigtHelper
<
dim
>::
size
;
FEEngine
&
interface_engine
=
emodel
.
getFEEngine
(
"EmbeddedInterfaceFEEngine"
);
const
auto
&
elem_filter
=
this
->
element_filter
(
interface_type
,
ghost_type
);
auto
nodes_per_background_e
=
Mesh
::
getNbNodesPerElement
(
background_type
);
auto
nb_quadrature_points
=
interface_engine
.
getNbIntegrationPoints
(
interface_type
,
ghost_type
);
auto
nb_element
=
elem_filter
.
size
();
auto
back_dof
=
dim
*
nodes_per_background_e
;
const
auto
&
shapesd
=
(
*
shape_derivatives
(
interface_type
,
ghost_type
))(
background_type
,
ghost_type
);
Array
<
Real
>
integrant
(
nb_quadrature_points
*
nb_element
,
back_dof
,
"integrant"
);
Matrix
<
Real
>
Bvoigt
(
voigt_size
,
back_dof
);
for
(
auto
&&
[
BtCt_sigma
,
C
,
sigma
,
B
]
:
zip
(
make_view
(
integrant
,
back_dof
),
make_view
(
directing_cosines
(
interface_type
,
ghost_type
),
voigt_size
),
this
->
stress
(
interface_type
,
ghost_type
),
make_view
(
shapesd
,
dim
,
nodes_per_background_e
)))
{
VoigtHelper
<
dim
>::
transferBMatrixToSymVoigtBMatrix
(
B
,
Bvoigt
,
nodes_per_background_e
);
BtCt_sigma
=
Bvoigt
.
transpose
()
*
C
*
sigma
*
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
<
Idx
>
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
,
Int
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
computeDirectingCosines
(
ElementType
type
,
GhostType
ghost_type
)
{
AKANTU_DEBUG_IN
();
Mesh
&
interface_mesh
=
emodel
.
getInterfaceMesh
();
const
auto
nb_nodes_per_element
=
Mesh
::
getNbNodesPerElement
(
type
);
const
auto
steel_dof
=
dim
*
nb_nodes_per_element
;
constexpr
auto
voigt_size
=
VoigtHelper
<
dim
>::
size
;
const
auto
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
));
for
(
auto
&&
data
:
zip
(
repeat_n
(
make_view
(
node_coordinates
,
dim
,
nb_nodes_per_element
),
nb_quad_points
),
make_view
<
1
,
voigt_size
>
(
directing_cosines
(
type
,
ghost_type
))))
{
computeDirectingCosinesOnQuad
(
std
::
get
<
0
>
(
data
),
std
::
get
<
1
>
(
data
));
}
// Mauro: the directing_cosines internal is defined on the quadrature points
// of each element
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
class
Mat
,
Int
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
assembleStiffnessMatrix
(
ElementType
type
,
GhostType
ghost_type
)
{
AKANTU_DEBUG_IN
();
Mesh
&
mesh
=
emodel
.
getMesh
();
for
(
auto
&&
mesh_type
:
mesh
.
elementTypes
(
dim
,
ghost_type
))
{
assembleStiffnessMatrixInterface
(
type
,
mesh_type
,
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
,
Int
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
assembleStiffnessMatrixInterface
(
ElementType
interface_type
,
ElementType
background_type
,
GhostType
ghost_type
)
{
AKANTU_DEBUG_IN
();
constexpr
Int
voigt_size
=
VoigtHelper
<
dim
>::
size
;
FEEngine
&
interface_engine
=
emodel
.
getFEEngine
(
"EmbeddedInterfaceFEEngine"
);
Array
<
Idx
>
&
elem_filter
=
this
->
element_filter
(
interface_type
,
ghost_type
);
Array
<
Real
>
&
grad_u
=
gradu_embedded
(
interface_type
,
ghost_type
);
Int
nb_element
=
elem_filter
.
size
();
Int
nodes_per_background_e
=
Mesh
::
getNbNodesPerElement
(
background_type
);
Int
nb_quadrature_points
=
interface_engine
.
getNbIntegrationPoints
(
interface_type
,
ghost_type
);
Int
back_dof
=
dim
*
nodes_per_background_e
;
Int
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
);
for
(
auto
&&
[
D
,
C
,
B
,
BtCtDCB
]
:
zip
(
tangent_moduli
,
make_view
(
directing_cosines
(
interface_type
,
ghost_type
),
1
,
voigt_size
),
make_view
(
shapesd
,
dim
,
nodes_per_background_e
),
make_view
(
integrant
,
integrant_size
,
integrant_size
)))
{
VoigtHelper
<
dim
>::
transferBMatrixToSymVoigtBMatrix
(
B
,
Bvoigt
,
nodes_per_background_e
);
BtCtDCB
=
Bvoigt
.
transpose
()
*
C
.
transpose
()
*
D
*
C
*
Bvoigt
*
area
;
}
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
,
Int
dim
>
Real
MaterialReinforcement
<
Mat
,
dim
>::
getEnergy
(
const
std
::
string
&
id
)
{
AKANTU_DEBUG_IN
();
if
(
id
==
"potential"
)
{
Real
epot
=
0.
;
this
->
computePotentialEnergyByElements
();
for
(
auto
&&
type
:
this
->
element_filter
.
elementTypes
(
this
->
spatial_dimension
))
{
FEEngine
&
interface_engine
=
emodel
.
getFEEngine
(
"EmbeddedInterfaceFEEngine"
);
epot
+=
interface_engine
.
integrate
(
this
->
potential_energy
(
type
,
_not_ghost
),
type
,
_not_ghost
,
this
->
element_filter
(
type
,
_not_ghost
));
epot
*=
area
;
}
return
epot
;
}
AKANTU_DEBUG_OUT
();
return
0
;
}
/* -------------------------------------------------------------------------- */
template
<
class
Mat
,
Int
dim
>
void
MaterialReinforcement
<
Mat
,
dim
>::
computeGradU
(
ElementType
interface_type
,
GhostType
ghost_type
)
{
// Looping over background types
for
(
auto
&&
bg_type
:
background_filter
(
interface_type
,
ghost_type
)
->
elementTypes
(
dim
))
{
const
Int
nodes_per_background_e
=
Mesh
::
getNbNodesPerElement
(
bg_type
);
constexpr
Int
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
,
dim
,
dim
>
concrete_du
(
dim
,
dim
);
Matrix
<
Real
,
dim
,
dim
>
epsilon
(
dim
,
dim
);
Vector
<
Real
,
voigt_size
>
evoigt
(
voigt_size
);
for
(
auto
&&
[
u
,
B
,
du
,
C
]
:
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
)))
{
concrete_du
=
u
*
B
.
transpose
();
auto
epsilon
=
0.5
*
(
concrete_du
+
concrete_du
.
transpose
());
VoigtHelper
<
dim
>::
matrixToVoigtWithFactors
(
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
,
Int
dim
>
template
<
class
Derived1
,
class
Derived2
>
inline
void
MaterialReinforcement
<
Mat
,
dim
>::
computeDirectingCosinesOnQuad
(
const
Eigen
::
MatrixBase
<
Derived1
>
&
nodes
,
Eigen
::
MatrixBase
<
Derived2
>
&
cosines
)
{
AKANTU_DEBUG_IN
();
AKANTU_DEBUG_ASSERT
(
nodes
.
cols
()
==
2
,
"Higher order reinforcement elements not implemented"
);
Vector
<
Real
>
delta
=
nodes
(
1
)
-
nodes
(
0
);
Real
sq_length
=
delta
.
dot
(
delta
);
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
();
}
/* -------------------------------------------------------------------------- */
}
// namespace akantu
Event Timeline
Log In to Comment