Homec4science
Diffusion akantu (bugfixes/cohesive_element_with_pbc) Browse

rAKA/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_tmpl.hh
50574cea5610bugfixes/cohesive_eleme…

Last Change
Raw File

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
c4science · Help