Homec4science
Diffusion akantu (bugfixes/60-hybrid-phase-field-formulation-unstable-in-dynamics) Browse

rAKA/src/fe_engine/shape_structural_inline_impl.hh
5acfa6ca970dbugfixes/60-hybrid…

Last Change
Raw File

shape_structural_inline_impl.hh
View Options

/**
* @file shape_structural_inline_impl.hh
*
* @author Fabian Barras <fabian.barras@epfl.ch>
* @author Lucas Frerot <lucas.frerot@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Wed Oct 11 2017
* @date last modification: Wed Feb 21 2018
*
* @brief ShapeStructural inline implementation
*
*
* Copyright (©) 2016-2018 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 "mesh_iterators.hh"
#include "shape_structural.hh"
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_SHAPE_STRUCTURAL_INLINE_IMPL_HH_
#define AKANTU_SHAPE_STRUCTURAL_INLINE_IMPL_HH_
namespace akantu {
namespace {
/// Extract nodal coordinates per elements
template <ElementType type>
std::unique_ptr<Array<Real>> getNodesPerElement(const Mesh & mesh,
const Array<Real> & nodes,
GhostType ghost_type) {
const auto dim = ElementClass<type>::getSpatialDimension();
const auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
auto nodes_per_element =
std::make_unique<Array<Real>>(0, dim * nb_nodes_per_element);
FEEngine::extractNodalToElementField(mesh, nodes, *nodes_per_element, type,
ghost_type);
return nodes_per_element;
}
} // namespace
template <ElementKind kind>
inline void ShapeStructural<kind>::initShapeFunctions(
const Array<Real> & /* unused */, const Matrix<Real> & /* unused */,
ElementType /* unused */, GhostType /* unused */) {
AKANTU_TO_IMPLEMENT();
}
/* -------------------------------------------------------------------------- */
#define INIT_SHAPE_FUNCTIONS(type) \
setIntegrationPointsByType<type>(integration_points, ghost_type); \
precomputeRotationMatrices<type>(nodes, ghost_type); \
precomputeShapesOnIntegrationPoints<type>(nodes, ghost_type); \
precomputeShapeDerivativesOnIntegrationPoints<type>(nodes, ghost_type);
template <>
inline void ShapeStructural<_ek_structural>::initShapeFunctions(
const Array<Real> & nodes, const Matrix<Real> & integration_points,
ElementType type, GhostType ghost_type) {
AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(INIT_SHAPE_FUNCTIONS);
}
#undef INIT_SHAPE_FUNCTIONS
/* -------------------------------------------------------------------------- */
template <ElementKind kind>
template <ElementType type>
void ShapeStructural<kind>::computeShapesOnIntegrationPointsInternal(
const Array<Real> & nodes, const Matrix<Real> & integration_points,
Array<Real> & shapes, GhostType ghost_type,
const Array<UInt> & filter_elements, bool mass) const {
auto nb_points = integration_points.cols();
auto nb_element = mesh.getConnectivity(type, ghost_type).size();
auto nb_nodes_per_element = ElementClass<type>::getNbNodesPerElement();
shapes.resize(nb_element * nb_points);
auto nb_dofs = ElementClass<type>::getNbDegreeOfFreedom();
auto nb_rows = nb_dofs;
if (mass) {
nb_rows = ElementClass<type>::getNbStressComponents();
}
#if !defined(AKANTU_NDEBUG)
UInt size_of_shapes = nb_rows * nb_dofs * nb_nodes_per_element;
AKANTU_DEBUG_ASSERT(shapes.getNbComponent() == size_of_shapes,
"The shapes array does not have the correct "
<< "number of component");
#endif
auto shapes_it = shapes.begin_reinterpret(
nb_rows, ElementClass<type>::getNbNodesPerInterpolationElement() * nb_dofs,
nb_points, nb_element);
auto shapes_begin = shapes_it;
if (filter_elements != empty_filter) {
nb_element = filter_elements.size();
}
auto nodes_per_element = getNodesPerElement<type>(mesh, nodes, ghost_type);
auto nodes_it = nodes_per_element->begin(mesh.getSpatialDimension(),
Mesh::getNbNodesPerElement(type));
auto nodes_begin = nodes_it;
auto rot_matrix_it =
make_view(rotation_matrices(type, ghost_type), nb_dofs, nb_dofs).begin();
auto rot_matrix_begin = rot_matrix_it;
for (UInt elem = 0; elem < nb_element; ++elem) {
if (filter_elements != empty_filter) {
shapes_it = shapes_begin + filter_elements(elem);
nodes_it = nodes_begin + filter_elements(elem);
rot_matrix_it = rot_matrix_begin + filter_elements(elem);
}
Tensor3<Real> & N = *shapes_it;
auto & real_coord = *nodes_it;
auto & RDOFs = *rot_matrix_it;
Matrix<Real> T(N.size(1), N.size(1), 0);
for (UInt i = 0; i < nb_nodes_per_element; ++i) {
T.block(RDOFs, i * RDOFs.rows(), i * RDOFs.rows());
}
if (not mass) {
ElementClass<type>::computeShapes(integration_points, real_coord, T, N);
} else {
ElementClass<type>::computeShapesMass(integration_points, real_coord, T,
N);
}
if (filter_elements == empty_filter) {
++shapes_it;
++nodes_it;
}
}
}
/* -------------------------------------------------------------------------- */
template <ElementKind kind>
template <ElementType type>
void ShapeStructural<kind>::precomputeRotationMatrices(
const Array<Real> & nodes, GhostType ghost_type) {
AKANTU_DEBUG_IN();
const auto spatial_dimension = mesh.getSpatialDimension();
const auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
const auto nb_element = mesh.getNbElement(type, ghost_type);
const auto nb_dof = ElementClass<type>::getNbDegreeOfFreedom();
if (not this->rotation_matrices.exists(type, ghost_type)) {
this->rotation_matrices.alloc(0, nb_dof * nb_dof, type, ghost_type);
}
auto & rot_matrices = this->rotation_matrices(type, ghost_type);
rot_matrices.resize(nb_element);
Array<Real> x_el(0, spatial_dimension * nb_nodes_per_element);
FEEngine::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type);
bool has_extra_normal = mesh.hasData<Real>("extra_normal", type, ghost_type);
Array<Real>::const_vector_iterator extra_normal;
if (has_extra_normal) {
extra_normal = mesh.getData<Real>("extra_normal", type, ghost_type)
.begin(spatial_dimension);
}
for (auto && tuple :
zip(make_view(x_el, spatial_dimension, nb_nodes_per_element),
make_view(rot_matrices, nb_dof, nb_dof))) {
// compute shape derivatives
auto & X = std::get<0>(tuple);
auto & R = std::get<1>(tuple);
if (has_extra_normal) {
ElementClass<type>::computeRotationMatrix(R, X, *extra_normal);
++extra_normal;
} else {
ElementClass<type>::computeRotationMatrix(
R, X, Vector<Real>(spatial_dimension));
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <ElementKind kind>
template <ElementType type>
void ShapeStructural<kind>::precomputeShapesOnIntegrationPoints(
const Array<Real> & nodes, GhostType ghost_type) {
AKANTU_DEBUG_IN();
const auto & natural_coords = integration_points(type, ghost_type);
auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
auto nb_points = integration_points(type, ghost_type).cols();
auto nb_element = mesh.getNbElement(type, ghost_type);
auto nb_dof = ElementClass<type>::getNbDegreeOfFreedom();
const auto dim = ElementClass<type>::getSpatialDimension();
const auto spatial_dimension = mesh.getSpatialDimension();
const auto natural_spatial_dimension =
ElementClass<type>::getNaturalSpaceDimension();
auto itp_type = FEEngine::getInterpolationType(type);
if (not shapes.exists(itp_type, ghost_type)) {
auto size_of_shapes = this->getShapeSize(type);
this->shapes.alloc(0, size_of_shapes, itp_type, ghost_type);
}
auto & rot_matrices = this->rotation_matrices(type, ghost_type);
auto & shapes_ = this->shapes(itp_type, ghost_type);
shapes_.resize(nb_element * nb_points);
auto nodes_per_element = getNodesPerElement<type>(mesh, nodes, ghost_type);
for (auto && tuple :
zip(make_view(shapes_, nb_dof, nb_dof * nb_nodes_per_element, nb_points),
make_view(*nodes_per_element, dim, nb_nodes_per_element),
make_view(rot_matrices, nb_dof, nb_dof))) {
auto & N = std::get<0>(tuple);
auto & X = std::get<1>(tuple);
auto & RDOFs = std::get<2>(tuple);
Matrix<Real> T(N.size(1), N.size(1), 0);
for (UInt i = 0; i < nb_nodes_per_element; ++i) {
T.block(RDOFs, i * RDOFs.rows(), i * RDOFs.rows());
}
auto R = RDOFs.block(0, 0, spatial_dimension, spatial_dimension);
// Rotate to local basis
auto x =
(R * X).block(0, 0, natural_spatial_dimension, nb_nodes_per_element);
ElementClass<type>::computeShapes(natural_coords, x, T, N);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <ElementKind kind>
template <ElementType type>
void ShapeStructural<kind>::precomputeShapeDerivativesOnIntegrationPoints(
const Array<Real> & nodes, GhostType ghost_type) {
AKANTU_DEBUG_IN();
const auto & natural_coords = integration_points(type, ghost_type);
const auto spatial_dimension = mesh.getSpatialDimension();
const auto natural_spatial_dimension =
ElementClass<type>::getNaturalSpaceDimension();
const auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
const auto nb_points = natural_coords.cols();
const auto nb_dof = ElementClass<type>::getNbDegreeOfFreedom();
const auto nb_element = mesh.getNbElement(type, ghost_type);
const auto nb_stress_components = ElementClass<type>::getNbStressComponents();
auto itp_type = FEEngine::getInterpolationType(type);
if (not this->shapes_derivatives.exists(itp_type, ghost_type)) {
auto size_of_shapesd = this->getShapeDerivativesSize(type);
this->shapes_derivatives.alloc(0, size_of_shapesd, itp_type, ghost_type);
}
auto & rot_matrices = this->rotation_matrices(type, ghost_type);
Array<Real> x_el(0, spatial_dimension * nb_nodes_per_element);
FEEngine::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type);
auto & shapesd = this->shapes_derivatives(itp_type, ghost_type);
shapesd.resize(nb_element * nb_points);
for (auto && tuple :
zip(make_view(x_el, spatial_dimension, nb_nodes_per_element),
make_view(shapesd, nb_stress_components,
nb_nodes_per_element * nb_dof, nb_points),
make_view(rot_matrices, nb_dof, nb_dof))) {
// compute shape derivatives
auto & X = std::get<0>(tuple);
auto & B = std::get<1>(tuple);
auto & RDOFs = std::get<2>(tuple);
Tensor3<Real> dnds(natural_spatial_dimension,
ElementClass<type>::interpolation_property::dnds_columns,
B.size(2));
ElementClass<type>::computeDNDS(natural_coords, X, dnds);
Tensor3<Real> J(natural_spatial_dimension, natural_spatial_dimension,
natural_coords.cols());
// Computing the coordinates of the element in the natural space
auto R = RDOFs.block(0, 0, spatial_dimension, spatial_dimension);
Matrix<Real> T(B.size(1), B.size(1), 0);
for (UInt i = 0; i < nb_nodes_per_element; ++i) {
T.block(RDOFs, i * RDOFs.rows(), i * RDOFs.rows());
}
// Rotate to local basis
auto x =
(R * X).block(0, 0, natural_spatial_dimension, nb_nodes_per_element);
ElementClass<type>::computeJMat(natural_coords, x, J);
ElementClass<type>::computeShapeDerivatives(J, dnds, T, B);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <ElementKind kind>
template <ElementType type>
void ShapeStructural<kind>::interpolateOnIntegrationPoints(
const Array<Real> & in_u, Array<Real> & out_uq, UInt nb_dof,
GhostType ghost_type, const Array<UInt> & filter_elements) const {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT(out_uq.getNbComponent() == nb_dof,
"The output array shape is not correct");
auto itp_type = FEEngine::getInterpolationType(type);
const auto & shapes_ = shapes(itp_type, ghost_type);
auto nb_element = mesh.getNbElement(type, ghost_type);
auto nb_nodes_per_element = ElementClass<type>::getNbNodesPerElement();
auto nb_quad_points_per_element = integration_points(type, ghost_type).cols();
Array<Real> u_el(0, nb_nodes_per_element * nb_dof);
FEEngine::extractNodalToElementField(mesh, in_u, u_el, type, ghost_type,
filter_elements);
auto nb_quad_points = nb_quad_points_per_element * u_el.size();
out_uq.resize(nb_quad_points);
auto out_it = out_uq.begin_reinterpret(nb_dof, 1, nb_quad_points_per_element,
u_el.size());
auto shapes_it =
shapes_.begin_reinterpret(nb_dof, nb_dof * nb_nodes_per_element,
nb_quad_points_per_element, nb_element);
auto u_it = u_el.begin_reinterpret(nb_dof * nb_nodes_per_element, 1,
nb_quad_points_per_element, u_el.size());
for_each_element(nb_element, filter_elements, [&](auto && el) {
auto & uq = *out_it;
const auto & u = *u_it;
auto N = Tensor3<Real>(shapes_it[el]);
for (auto && q : arange(uq.size(2))) {
auto uq_q = Matrix<Real>(uq(q));
auto u_q = Matrix<Real>(u(q));
auto N_q = Matrix<Real>(N(q));
uq_q.mul<false, false>(N_q, u_q);
}
++out_it;
++u_it;
});
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <ElementKind kind>
template <ElementType type>
void ShapeStructural<kind>::gradientOnIntegrationPoints(
const Array<Real> & in_u, Array<Real> & out_nablauq, UInt nb_dof,
GhostType ghost_type, const Array<UInt> & filter_elements) const {
AKANTU_DEBUG_IN();
auto itp_type = FEEngine::getInterpolationType(type);
const auto & shapesd = shapes_derivatives(itp_type, ghost_type);
auto nb_element = mesh.getNbElement(type, ghost_type);
auto element_dimension = ElementClass<type>::getSpatialDimension();
auto nb_quad_points_per_element = integration_points(type, ghost_type).cols();
auto nb_nodes_per_element = ElementClass<type>::getNbNodesPerElement();
Array<Real> u_el(0, nb_nodes_per_element * nb_dof);
FEEngine::extractNodalToElementField(mesh, in_u, u_el, type, ghost_type,
filter_elements);
auto nb_quad_points = nb_quad_points_per_element * u_el.size();
out_nablauq.resize(nb_quad_points);
auto out_it = out_nablauq.begin_reinterpret(
element_dimension, 1, nb_quad_points_per_element, u_el.size());
auto shapesd_it = shapesd.begin_reinterpret(
element_dimension, nb_dof * nb_nodes_per_element,
nb_quad_points_per_element, nb_element);
auto u_it = u_el.begin_reinterpret(nb_dof * nb_nodes_per_element, 1,
nb_quad_points_per_element, u_el.size());
for_each_element(nb_element, filter_elements, [&](auto && el) {
auto & nablau = *out_it;
const auto & u = *u_it;
auto B = Tensor3<Real>(shapesd_it[el]);
for (auto && q : arange(nablau.size(2))) {
auto nablau_q = Matrix<Real>(nablau(q));
auto u_q = Matrix<Real>(u(q));
auto B_q = Matrix<Real>(B(q));
nablau_q.mul<false, false>(B_q, u_q);
}
++out_it;
++u_it;
});
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <>
template <ElementType type>
void ShapeStructural<_ek_structural>::computeBtD(
const Array<Real> & Ds, Array<Real> & BtDs, GhostType ghost_type,
const Array<UInt> & filter_elements) const {
auto itp_type = ElementClassProperty<type>::interpolation_type;
auto nb_stress = ElementClass<type>::getNbStressComponents();
auto nb_dof_per_element = ElementClass<type>::getNbDegreeOfFreedom() *
mesh.getNbNodesPerElement(type);
const auto & shapes_derivatives =
this->shapes_derivatives(itp_type, ghost_type);
Array<Real> shapes_derivatives_filtered(0,
shapes_derivatives.getNbComponent());
auto && view = make_view(shapes_derivatives, nb_stress, nb_dof_per_element);
auto B_it = view.begin();
auto B_end = view.end();
if (filter_elements != empty_filter) {
FEEngine::filterElementalData(this->mesh, shapes_derivatives,
shapes_derivatives_filtered, type, ghost_type,
filter_elements);
auto && view =
make_view(shapes_derivatives_filtered, nb_stress, nb_dof_per_element);
B_it = view.begin();
B_end = view.end();
}
for (auto && values : zip(range(B_it, B_end), make_view(Ds, nb_stress),
make_view(BtDs, BtDs.getNbComponent()))) {
const auto & B = std::get<0>(values);
const auto & D = std::get<1>(values);
auto & Bt_D = std::get<2>(values);
Bt_D.template mul<true>(B, D);
}
}
/* -------------------------------------------------------------------------- */
template <>
template <ElementType type>
void ShapeStructural<_ek_structural>::computeNtb(
const Array<Real> & bs, Array<Real> & Ntbs, GhostType ghost_type,
const Array<UInt> & filter_elements) const {
auto itp_type = ElementClassProperty<type>::interpolation_type;
auto nb_dof = ElementClass<type>::getNbDegreeOfFreedom();
auto nb_nodes_per_element = mesh.getNbNodesPerElement(type);
const auto & shapes = this->shapes(itp_type, ghost_type);
Array<Real> shapes_filtered(0, shapes.getNbComponent());
auto && view = make_view(shapes, nb_dof, nb_dof * nb_nodes_per_element);
auto N_it = view.begin();
auto N_end = view.end();
if (filter_elements != empty_filter) {
FEEngine::filterElementalData(this->mesh, shapes, shapes_filtered, type,
ghost_type, filter_elements);
auto && view =
make_view(shapes_filtered, nb_dof, nb_dof * nb_nodes_per_element);
N_it = view.begin();
N_end = view.end();
}
for (auto && values : zip(range(N_it, N_end), make_view(bs, nb_dof),
make_view(Ntbs, nb_dof * nb_nodes_per_element))) {
const auto & N = std::get<0>(values);
const auto & b = std::get<1>(values);
auto & Nt_b = std::get<2>(values);
Nt_b.template mul<true>(N, b);
}
}
} // namespace akantu
#endif /* AKANTU_SHAPE_STRUCTURAL_INLINE_IMPL_HH_ */
c4science · Help