Homec4science
Diffusion akantu (bugfixes/39-the-python-example-for-phasefield-are-broken-due-to-a-change-to-setintegrationscheme) Browse

rAKA/src/fe_engine/gauss_integration_tmpl.hh
d213b95f146fbugfixes/39-the-python…

Last Change
Raw File

gauss_integration_tmpl.hh
View Options

/**
* @file gauss_integration_tmpl.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Tue May 10 2016
* @date last modification: Tue Sep 29 2020
*
* @brief implementation of the gauss integration helpers
*
*
* @section LICENSE
*
* Copyright (©) 2016-2021 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/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_GAUSS_INTEGRATION_TMPL_HH_
#define AKANTU_GAUSS_INTEGRATION_TMPL_HH_
namespace akantu {
/* -------------------------------------------------------------------------- */
/* GaussIntegrationElement */
/* -------------------------------------------------------------------------- */
namespace _aka_gauss_helpers {
template <GaussIntegrationType type, UInt n> struct GaussIntegrationNbPoints {
static const UInt nb_points = 0;
};
#if !defined(DOXYGEN)
template <UInt n> struct GaussIntegrationNbPoints<_git_not_defined, n> {
static const UInt nb_points = 0;
};
template <UInt n> struct GaussIntegrationNbPoints<_git_point, n> {
static const UInt nb_points = 1;
};
template <UInt n> struct GaussIntegrationNbPoints<_git_segment, n> {
static const UInt nb_points = (n + 1) / 2 + (bool((n + 1) % 2) ? 1 : 0);
};
#define DECLARE_GAUSS_NB_POINTS(type, order, points) \
template <> struct GaussIntegrationNbPoints<type, order> { \
static const UInt nb_points = points; \
}
#define DECLARE_GAUSS_NB_POINTS_PENT(type, order, xo, yo) \
template <> struct GaussIntegrationNbPoints<type, order> { \
static const UInt x_order = xo; \
static const UInt yz_order = yo; \
static const UInt nb_points = 1; \
}
DECLARE_GAUSS_NB_POINTS(_git_triangle, 1, 1);
DECLARE_GAUSS_NB_POINTS(_git_triangle, 2, 3);
DECLARE_GAUSS_NB_POINTS(_git_triangle, 3, 4);
DECLARE_GAUSS_NB_POINTS(_git_triangle, 4, 6);
DECLARE_GAUSS_NB_POINTS(_git_triangle, 5, 7);
DECLARE_GAUSS_NB_POINTS(_git_tetrahedron, 1, 1);
DECLARE_GAUSS_NB_POINTS(_git_tetrahedron, 2, 4);
DECLARE_GAUSS_NB_POINTS(_git_tetrahedron, 3, 5);
DECLARE_GAUSS_NB_POINTS(_git_tetrahedron, 4, 15);
DECLARE_GAUSS_NB_POINTS(_git_tetrahedron, 5, 15);
DECLARE_GAUSS_NB_POINTS_PENT(_git_pentahedron, 1, 3,
2); // order 3 in x, order 2 in y and z
DECLARE_GAUSS_NB_POINTS_PENT(_git_pentahedron, 2, 3,
2); // order 3 in x, order 2 in y and z
DECLARE_GAUSS_NB_POINTS_PENT(_git_pentahedron, 3, 3,
3); // order 3 in x, order 3 in y and z
DECLARE_GAUSS_NB_POINTS_PENT(_git_pentahedron, 4, 5,
5); // order 5 in x, order 5 in y and z
DECLARE_GAUSS_NB_POINTS_PENT(_git_pentahedron, 5, 5,
5); // order 5 in x, order 5 in y and z
template <GaussIntegrationType type, UInt n, UInt on = n,
bool end_recurse = false>
struct GaussIntegrationNbPointsHelper {
static const UInt pnp = GaussIntegrationNbPoints<type, n>::nb_points;
static const UInt order = n;
static const UInt nb_points = pnp;
};
template <GaussIntegrationType type, UInt n, UInt on>
struct GaussIntegrationNbPointsHelper<type, n, on, true> {
static const UInt nb_points = 0;
};
#endif
/* ------------------------------------------------------------------------ */
/* Generic helper */
/* ------------------------------------------------------------------------ */
template <GaussIntegrationType type, UInt dimension, UInt n>
struct GaussIntegrationTypeDataHelper {
using git_np = GaussIntegrationNbPoints<type, n>;
using git_data = GaussIntegrationTypeData<type, git_np::nb_points>;
static UInt getNbQuadraturePoints() { return git_np::nb_points; }
static Matrix<Real> getQuadraturePoints() {
return Matrix<Real>(git_data::quad_positions, dimension,
git_np::nb_points);
}
static Vector<Real> getWeights() {
return Vector<Real>(git_data::quad_weights, git_np::nb_points);
}
};
#if !defined(DOXYGEN)
/* ------------------------------------------------------------------------ */
/* helper for _segment _quadrangle _hexahedron */
/* ------------------------------------------------------------------------ */
template <UInt dimension, UInt dp>
struct GaussIntegrationTypeDataHelper<_git_segment, dimension, dp> {
using git_np = GaussIntegrationNbPoints<_git_segment, dp>;
using git_data = GaussIntegrationTypeData<_git_segment, git_np::nb_points>;
static UInt getNbQuadraturePoints() {
return Math::pow<dimension>(git_np::nb_points);
}
static Matrix<Real> getQuadraturePoints() {
UInt tot_nquad = getNbQuadraturePoints();
UInt nquad = git_np::nb_points;
Matrix<Real> quads(dimension, tot_nquad);
Vector<Real> pos(git_data::quad_positions, nquad);
UInt offset = 1;
for (UInt d = 0; d < dimension; ++d) {
for (UInt n = 0, q = 0; n < tot_nquad; ++n, q += offset) {
UInt rq = q % tot_nquad + q / tot_nquad;
quads(d, rq) = pos(n % nquad);
}
offset *= nquad;
}
return quads;
}
static Vector<Real> getWeights() {
UInt tot_nquad = getNbQuadraturePoints();
UInt nquad = git_np::nb_points;
Vector<Real> quads_weights(tot_nquad, 1.);
Vector<Real> weights(git_data::quad_weights, nquad);
UInt offset = 1;
for (UInt d = 0; d < dimension; ++d) {
for (UInt n = 0, q = 0; n < tot_nquad; ++n, q += offset) {
UInt rq = q % tot_nquad + q / tot_nquad;
quads_weights(rq) *= weights(n % nquad);
}
offset *= nquad;
}
return quads_weights;
}
};
/* ------------------------------------------------------------------------ */
/* helper for _pentahedron */
/* ------------------------------------------------------------------------ */
template <UInt dimension, UInt dp>
struct GaussIntegrationTypeDataHelper<_git_pentahedron, dimension, dp> {
using git_info = GaussIntegrationNbPoints<_git_pentahedron, dp>;
using git_np_seg =
GaussIntegrationNbPoints<_git_segment, git_info::x_order>;
using git_np_tri =
GaussIntegrationNbPoints<_git_triangle, git_info::yz_order>;
using git_data_seg =
GaussIntegrationTypeData<_git_segment, git_np_seg::nb_points>;
using git_data_tri =
GaussIntegrationTypeData<_git_triangle, git_np_tri::nb_points>;
static UInt getNbQuadraturePoints() {
return git_np_seg::nb_points * git_np_tri::nb_points;
}
static Matrix<Real> getQuadraturePoints() {
UInt tot_nquad = getNbQuadraturePoints();
UInt nquad_seg = git_np_seg::nb_points;
UInt nquad_tri = git_np_tri::nb_points;
Matrix<Real> quads(dimension, tot_nquad);
Matrix<Real> pos_seg_w(git_data_seg::quad_positions, 1, nquad_seg);
Matrix<Real> pos_tri_w(git_data_tri::quad_positions, 2, nquad_tri);
for (UInt ns = 0, q = 0; ns < nquad_seg; ++ns) {
Vector<Real> pos_seg = pos_seg_w(ns);
for (UInt nt = 0; nt < nquad_tri; ++nt, ++q) {
Vector<Real> pos_tri = pos_tri_w(nt);
Vector<Real> quad = quads(q);
quad(_x) = pos_seg(_x);
quad(_y) = pos_tri(_x);
quad(_z) = pos_tri(_y);
}
}
return quads;
}
static Vector<Real> getWeights() {
UInt tot_nquad = getNbQuadraturePoints();
UInt nquad_seg = git_np_seg::nb_points;
UInt nquad_tri = git_np_tri::nb_points;
Vector<Real> quads_weights(tot_nquad);
Vector<Real> weight_seg(git_data_seg::quad_weights, nquad_seg);
Vector<Real> weight_tri(git_data_tri::quad_weights, nquad_tri);
for (UInt ns = 0, q = 0; ns < nquad_seg; ++ns) {
for (UInt nt = 0; nt < nquad_tri; ++nt, ++q) {
quads_weights(q) = weight_seg(ns) * weight_tri(nt);
}
}
return quads_weights;
}
};
#endif
} // namespace _aka_gauss_helpers
template <ElementType element_type, UInt n>
Matrix<Real> GaussIntegrationElement<element_type, n>::getQuadraturePoints() {
const InterpolationType itp_type =
ElementClassProperty<element_type>::interpolation_type;
using interpolation_property = InterpolationProperty<itp_type>;
using data_helper = _aka_gauss_helpers::GaussIntegrationTypeDataHelper<
ElementClassProperty<element_type>::gauss_integration_type,
interpolation_property::natural_space_dimension, n>;
Matrix<Real> tmp(data_helper::getQuadraturePoints());
return tmp;
}
/* -------------------------------------------------------------------------- */
template <ElementType element_type, UInt n>
Vector<Real> GaussIntegrationElement<element_type, n>::getWeights() {
const InterpolationType itp_type =
ElementClassProperty<element_type>::interpolation_type;
using interpolation_property = InterpolationProperty<itp_type>;
using data_helper = _aka_gauss_helpers::GaussIntegrationTypeDataHelper<
ElementClassProperty<element_type>::gauss_integration_type,
interpolation_property::natural_space_dimension, n>;
Vector<Real> tmp(data_helper::getWeights());
return tmp;
}
/* -------------------------------------------------------------------------- */
template <ElementType element_type, UInt n>
UInt GaussIntegrationElement<element_type, n>::getNbQuadraturePoints() {
const InterpolationType itp_type =
ElementClassProperty<element_type>::interpolation_type;
using interpolation_property = InterpolationProperty<itp_type>;
using data_helper = _aka_gauss_helpers::GaussIntegrationTypeDataHelper<
ElementClassProperty<element_type>::gauss_integration_type,
interpolation_property::natural_space_dimension, n>;
return data_helper::getNbQuadraturePoints();
}
} // namespace akantu
#endif /* AKANTU_GAUSS_INTEGRATION_TMPL_HH_ */
c4science · Help