diff --git a/src/fe_engine/fe_engine_template_tmpl_field.hh b/src/fe_engine/fe_engine_template_tmpl_field.hh
index 225c145ec..5b8641797 100644
--- a/src/fe_engine/fe_engine_template_tmpl_field.hh
+++ b/src/fe_engine/fe_engine_template_tmpl_field.hh
@@ -1,428 +1,426 @@
/**
* @file fe_engine_template_tmpl_field.hh
*
* @author Nicolas Richart
*
* @date creation Tue Jul 25 2017
*
* @brief implementation of the assemble field s functions
*
* @section LICENSE
*
* Copyright (©) 2010-2011 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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "fe_engine_template.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_FE_ENGINE_TEMPLATE_TMPL_FIELD_HH__
#define __AKANTU_FE_ENGINE_TEMPLATE_TMPL_FIELD_HH__
namespace akantu {
/* -------------------------------------------------------------------------- */
/* Matrix lumping functions */
/* -------------------------------------------------------------------------- */
namespace fe_engine {
namespace details {
namespace {
template
void fillField(const Functor & field_funct, Array & field,
UInt nb_element, UInt nb_integration_points,
const ElementType & type, const GhostType & ghost_type) {
UInt nb_degree_of_freedom = field.getNbComponent();
field.resize(nb_integration_points * nb_element);
auto field_it = field.begin_reinterpret(
nb_degree_of_freedom, nb_integration_points, nb_element);
Element el(type, 0, ghost_type);
for (; el.element < nb_element; ++el.element, ++field_it) {
field_funct(*field_it, el);
}
}
} // namespace
} // namespace details
} // namespace fe_engine
/**
* Helper class to be able to write a partial specialization on the element kind
*/
namespace fe_engine {
namespace details {
template struct AssembleLumpedTemplateHelper {};
#define ASSEMBLE_LUMPED(type) \
fem.template assembleFieldLumped(field_funct, lumped, dof_id, \
dof_manager, ghost_type)
#define AKANTU_SPECIALIZE_ASSEMBLE_HELPER(kind) \
template <> struct AssembleLumpedTemplateHelper { \
template class I, \
template class S, ElementKind k, class IOF, \
class Functor> \
static void call(const FEEngineTemplate & fem, \
const Functor & field_funct, const ID & lumped, \
const ID & dof_id, DOFManager & dof_manager, \
ElementType type, const GhostType & ghost_type) { \
AKANTU_BOOST_KIND_ELEMENT_SWITCH(ASSEMBLE_LUMPED, kind); \
} \
};
AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_ASSEMBLE_HELPER)
#undef AKANTU_SPECIALIZE_ASSEMBLE_HELPER
#undef ASSEMBLE_LUMPED
} // namespace details
} // namespace fe_engine
/* -------------------------------------------------------------------------- */
template class I, template class S,
ElementKind kind, class IOF>
template
void FEEngineTemplate::assembleFieldLumped(
const Functor & field_funct, const ID & matrix_id, const ID & dof_id,
DOFManager & dof_manager, ElementType type,
const GhostType & ghost_type) const {
AKANTU_DEBUG_IN();
fe_engine::details::AssembleLumpedTemplateHelper::call(
*this, field_funct, matrix_id, dof_id, dof_manager, type, ghost_type);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
template class I, template class S,
ElementKind kind, class IntegrationOrderFunctor>
template
void FEEngineTemplate::assembleFieldLumped(
const Functor & field_funct, const ID & matrix_id, const ID & dof_id,
DOFManager & dof_manager, const GhostType & ghost_type) const {
AKANTU_DEBUG_IN();
UInt nb_degree_of_freedom = dof_manager.getDOFs(dof_id).getNbComponent();
UInt nb_element = mesh.getNbElement(type, ghost_type);
UInt nb_integration_points = this->getNbIntegrationPoints(type);
Array field(0, nb_degree_of_freedom);
fe_engine::details::fillField(field_funct, field, nb_element,
nb_integration_points, type, ghost_type);
switch (type) {
case _triangle_6:
case _quadrangle_8:
case _tetrahedron_10:
case _hexahedron_20:
case _pentahedron_15:
this->template assembleLumpedDiagonalScaling(field, matrix_id, dof_id,
dof_manager, ghost_type);
+ break;
default:
this->template assembleLumpedRowSum(field, matrix_id, dof_id,
dof_manager, ghost_type);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* @f$ \tilde{M}_{i} = \sum_j M_{ij} = \sum_j \int \rho \varphi_i \varphi_j dV =
* \int \rho \varphi_i dV @f$
*/
template class I, template class S,
ElementKind kind, class IntegrationOrderFunctor>
template
void FEEngineTemplate::
assembleLumpedRowSum(const Array & field, const ID & matrix_id,
const ID & dof_id, DOFManager & dof_manager,
const GhostType & ghost_type) const {
AKANTU_DEBUG_IN();
UInt shapes_size = ElementClass::getShapeSize();
UInt nb_degree_of_freedom = field.getNbComponent();
Array * field_times_shapes =
new Array(0, shapes_size * nb_degree_of_freedom);
shape_functions.template fieldTimesShapes(field, *field_times_shapes,
ghost_type);
UInt nb_element = mesh.getNbElement(type, ghost_type);
Array * int_field_times_shapes = new Array(
nb_element, shapes_size * nb_degree_of_freedom, "inte_rho_x_shapes");
integrator.template integrate(
*field_times_shapes, *int_field_times_shapes,
nb_degree_of_freedom * shapes_size, ghost_type, empty_filter);
delete field_times_shapes;
dof_manager.assembleElementalArrayToLumpedMatrix(
dof_id, *int_field_times_shapes, matrix_id, type, ghost_type);
delete int_field_times_shapes;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* @f$ \tilde{M}_{i} = c * M_{ii} = \int_{V_e} \rho dV @f$
*/
template class I, template class S,
ElementKind kind, class IntegrationOrderFunctor>
template
void FEEngineTemplate::
assembleLumpedDiagonalScaling(const Array & field,
const ID & matrix_id, const ID & dof_id,
DOFManager & dof_manager,
const GhostType & ghost_type) const {
AKANTU_DEBUG_IN();
const ElementType & type_p1 = ElementClass::getP1ElementType();
UInt nb_nodes_per_element_p1 = Mesh::getNbNodesPerElement(type_p1);
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
UInt nb_degree_of_freedom = field.getNbComponent();
UInt nb_element = mesh.getNbElement(type, ghost_type);
Vector nodal_factor(nb_nodes_per_element);
#define ASSIGN_WEIGHT_TO_NODES(corner, mid) \
{ \
for (UInt n = 0; n < nb_nodes_per_element_p1; n++) \
nodal_factor(n) = corner; \
for (UInt n = nb_nodes_per_element_p1; n < nb_nodes_per_element; n++) \
nodal_factor(n) = mid; \
}
if (type == _triangle_6)
ASSIGN_WEIGHT_TO_NODES(1. / 12., 1. / 4.);
if (type == _tetrahedron_10)
ASSIGN_WEIGHT_TO_NODES(1. / 32., 7. / 48.);
if (type == _quadrangle_8)
ASSIGN_WEIGHT_TO_NODES(
3. / 76.,
16. / 76.); /** coeff. derived by scaling
* the diagonal terms of the corresponding
* consistent mass computed with 3x3 gauss points;
* coeff. are (1./36., 8./36.) for 2x2 gauss points */
if (type == _hexahedron_20)
ASSIGN_WEIGHT_TO_NODES(
7. / 248.,
16. / 248.); /** coeff. derived by scaling
* the diagonal terms of the corresponding
* consistent mass computed with 3x3x3 gauss points;
* coeff. are (1./40.,
* 1./15.) for 2x2x2 gauss points */
if (type == _pentahedron_15) {
// coefficients derived by scaling the diagonal terms of the corresponding
// consistent mass computed with 8 gauss points;
for (UInt n = 0; n < nb_nodes_per_element_p1; n++)
nodal_factor(n) = 51. / 2358.;
Real mid_triangle = 192. / 2358.;
Real mid_quadrangle = 300. / 2358.;
nodal_factor(6) = mid_triangle;
nodal_factor(7) = mid_triangle;
nodal_factor(8) = mid_triangle;
nodal_factor(9) = mid_quadrangle;
nodal_factor(10) = mid_quadrangle;
nodal_factor(11) = mid_quadrangle;
nodal_factor(12) = mid_triangle;
nodal_factor(13) = mid_triangle;
nodal_factor(14) = mid_triangle;
}
if (nb_element == 0) {
AKANTU_DEBUG_OUT();
return;
}
#undef ASSIGN_WEIGHT_TO_NODES
-
/// compute @f$ \int \rho dV = \rho V @f$ for each element
- Array * int_field =
- new Array(field.size(), nb_degree_of_freedom, "inte_rho_x");
+ auto int_field =
+ std::make_unique>(field.size(), nb_degree_of_freedom, "inte_rho_x");
integrator.template integrate(field, *int_field, nb_degree_of_freedom,
ghost_type, empty_filter);
/// distribute the mass of the element to the nodes
- Array * lumped_per_node = new Array(
+ auto lumped_per_node = std::make_unique>(
nb_element, nb_degree_of_freedom * nb_nodes_per_element, "mass_per_node");
- Array::const_vector_iterator int_field_it =
+ auto int_field_it =
int_field->begin(nb_degree_of_freedom);
- Array::matrix_iterator lumped_per_node_it =
+ auto lumped_per_node_it =
lumped_per_node->begin(nb_degree_of_freedom, nb_nodes_per_element);
for (UInt e = 0; e < nb_element; ++e) {
for (UInt n = 0; n < nb_nodes_per_element; ++n) {
Vector l = (*lumped_per_node_it)(n);
l = *int_field_it;
l *= nodal_factor(n);
}
++int_field_it;
++lumped_per_node_it;
}
- delete int_field;
dof_manager.assembleElementalArrayToLumpedMatrix(dof_id, *lumped_per_node,
matrix_id, type, ghost_type);
- delete lumped_per_node;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* Helper class to be able to write a partial specialization on the element kind
*/
namespace fe_engine {
namespace details {
template struct AssembleFieldMatrixHelper {};
#define ASSEMBLE_MATRIX(type) \
fem.template assembleFieldMatrix( \
field_funct, matrix_id, dof_id, dof_manager, ghost_type)
#define AKANTU_SPECIALIZE_ASSEMBLE_FIELD_MATRIX_HELPER(kind) \
template <> struct AssembleFieldMatrixHelper { \
template class I, \
template class S, ElementKind k, class IOF, \
class Functor> \
static void call(const FEEngineTemplate & fem, \
const Functor & field_funct, const ID & matrix_id, \
const ID & dof_id, DOFManager & dof_manager, \
ElementType type, const GhostType & ghost_type) { \
AKANTU_BOOST_KIND_ELEMENT_SWITCH(ASSEMBLE_MATRIX, kind); \
} \
};
AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_ASSEMBLE_FIELD_MATRIX_HELPER)
#undef AKANTU_SPECIALIZE_ASSEMBLE_FIELD_MATRIX_HELPER
#undef ASSEMBLE_MATRIX
} // namespace details
} // namespace fe_engine
/* -------------------------------------------------------------------------- */
template class I, template class S,
ElementKind kind, class IOF>
template
void FEEngineTemplate::assembleFieldMatrix(
const Functor & field_funct, const ID & matrix_id, const ID & dof_id,
DOFManager & dof_manager, ElementType type,
const GhostType & ghost_type) const {
AKANTU_DEBUG_IN();
fe_engine::details::AssembleFieldMatrixHelper::template call(
*this, field_funct, matrix_id, dof_id, dof_manager, type, ghost_type);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* @f$ \tilde{M}_{i} = \sum_j M_{ij} = \sum_j \int \rho \varphi_i \varphi_j dV =
* \int \rho \varphi_i dV @f$
*/
template class I, template class S,
ElementKind kind, class IntegrationOrderFunctor>
template
void FEEngineTemplate::assembleFieldMatrix(
const Functor & field_funct, const ID & matrix_id, const ID & dof_id,
DOFManager & dof_manager, const GhostType & ghost_type) const {
AKANTU_DEBUG_IN();
UInt shapes_size = ElementClass::getShapeSize();
UInt nb_degree_of_freedom = dof_manager.getDOFs(dof_id).getNbComponent();
UInt lmat_size = nb_degree_of_freedom * shapes_size;
UInt nb_element = mesh.getNbElement(type, ghost_type);
// \int N * N so degree 2 * degree of N
const UInt polynomial_degree =
2 * ElementClassProperty::polynomial_degree;
// getting the integration points
Matrix integration_points =
integrator.template getIntegrationPoints();
UInt nb_integration_points = integration_points.cols();
UInt vect_size = nb_integration_points * nb_element;
// getting the shapes on the integration points
Array shapes(0, shapes_size);
shape_functions.template computeShapesOnIntegrationPoints(
mesh.getNodes(), integration_points, shapes, ghost_type);
// Extending the shape functions
Array modified_shapes(vect_size, lmat_size * nb_degree_of_freedom, 0.);
Array local_mat(vect_size, lmat_size * lmat_size);
auto mshapes_it = modified_shapes.begin(nb_degree_of_freedom, lmat_size);
auto shapes_it = shapes.begin(shapes_size);
for (UInt q = 0; q < vect_size; ++q, ++mshapes_it, ++shapes_it) {
for (UInt d = 0; d < nb_degree_of_freedom; ++d) {
for (UInt s = 0; s < shapes_size; ++s) {
(*mshapes_it)(d, s * nb_degree_of_freedom + d) = (*shapes_it)(s);
}
}
}
// getting the value to assemble on the integration points
Array field(vect_size, nb_degree_of_freedom);
fe_engine::details::fillField(field_funct, field, nb_element,
nb_integration_points, type, ghost_type);
// computing \rho * N
mshapes_it = modified_shapes.begin(nb_degree_of_freedom, lmat_size);
auto lmat = local_mat.begin(lmat_size, lmat_size);
auto field_it =
field.begin_reinterpret(nb_degree_of_freedom, field.size());
for (UInt q = 0; q < vect_size; ++q, ++lmat, ++mshapes_it, ++field_it) {
const auto & rho = *field_it;
const auto & N = *mshapes_it;
auto & mat = *lmat;
Matrix Nt = N.transpose();
for (UInt d = 0; d < Nt.cols(); ++d) {
Nt(d) *= rho(d);
}
mat.mul(Nt, N);
}
// integrate the elemental values
Array int_field_times_shapes(nb_element, lmat_size * lmat_size,
"inte_rho_x_shapes");
this->integrator.template integrate(
local_mat, int_field_times_shapes, lmat_size * lmat_size, ghost_type);
// assemble the elemental values to the matrix
dof_manager.assembleElementalMatricesToMatrix(
matrix_id, dof_id, int_field_times_shapes, type, ghost_type);
AKANTU_DEBUG_OUT();
}
} // namespace akantu
#endif /* __AKANTU_FE_ENGINE_TEMPLATE_TMPL_FIELD_HH__ */
diff --git a/test/test_model/test_solid_mechanics_model/patch_tests/lumped_mass/test_lumped_mass.cc b/test/test_model/test_solid_mechanics_model/patch_tests/lumped_mass/test_lumped_mass.cc
index 77c4d531d..da518ef24 100644
--- a/test/test_model/test_solid_mechanics_model/patch_tests/lumped_mass/test_lumped_mass.cc
+++ b/test/test_model/test_solid_mechanics_model/patch_tests/lumped_mass/test_lumped_mass.cc
@@ -1,91 +1,90 @@
/**
* @file test_lumped_mass.cc
*
* @author Daniel Pino Muñoz
* @author Nicolas Richart
*
* @date creation: Thu Mar 27 2014
* @date last modification: Thu Oct 15 2015
*
* @brief test the lumping of the mass matrix
*
* @section LICENSE
*
* Copyright (©) 2014, 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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "solid_mechanics_model.hh"
using namespace akantu;
/* -------------------------------------------------------------------------- */
int main(int argc, char *argv[])
{
UInt spatial_dimension = ElementClass::getSpatialDimension();
const ElementType type = TYPE;
akantu::initialize("material.dat", argc, argv);
Mesh mesh(spatial_dimension);
std::stringstream filename; filename << TYPE << ".msh";
mesh.read(filename.str());
SolidMechanicsModel model(mesh);
/// model initialization
model.initFull();
model.assembleMassLumped();
Real rho = model.getMaterial(0).getRho();
FEEngine & fem = model.getFEEngine();
UInt nb_element = mesh.getNbElement(type);
UInt nb_quadrature_points = fem.getNbIntegrationPoints(type) * nb_element;
Array rho_on_quad(nb_quadrature_points, 1, rho, "rho_on_quad");
Real mass = fem.integrate(rho_on_quad, type);
Array::const_vector_iterator mass_it = model.getMass().begin(spatial_dimension);
Array::const_vector_iterator mass_end = model.getMass().end(spatial_dimension);
Vector sum(spatial_dimension, 0.);
-
for(; mass_it != mass_end; ++mass_it) {
sum += *mass_it;
}
std::cout << mass << std::endl << sum << std::endl;
if(!(std::abs((mass - sum[0])/mass) < 2e-15)) {
std::cerr << "total mass is not correct" << std::endl;
return EXIT_FAILURE;
}
for(UInt i = 1; i < spatial_dimension; ++i)
if(!(std::abs((sum[0] - sum[i])/mass) < 1e-15)) {
std::cerr << "total mass is not correct" << std::endl;
return EXIT_FAILURE;
}
finalize();
return EXIT_SUCCESS;
}