diff --git a/ProcessLib/Deformation/LinearBMatrix.h b/ProcessLib/Deformation/LinearBMatrix.h
index 1d122c9a14244f4df6cb8fb60f9eeb419907c1ed..139d5c0c85cd45c9e394682e1d291bbc2227c1b4 100644
--- a/ProcessLib/Deformation/LinearBMatrix.h
+++ b/ProcessLib/Deformation/LinearBMatrix.h
@@ -18,13 +18,14 @@ namespace LinearBMatrix
namespace detail
{
template <int NPOINTS, typename DNDX_Type, typename BMatrixType>
-void fillBMatrix2DCartesianPart(DNDX_Type const& dNdx, BMatrixType& b_matrix)
+void fillBMatrix2DCartesianPart(DNDX_Type const& dNdx, BMatrixType& B)
{
- for (int i = 0; i < NPOINTS; ++i) {
- b_matrix(1, NPOINTS + i) = dNdx(1, i);
- b_matrix(3, i) = dNdx(1, i) / std::sqrt(2);
- b_matrix(3, NPOINTS + i) = dNdx(0, i) / std::sqrt(2);
- b_matrix(0, i) = dNdx(0, i);
+ for (int i = 0; i < NPOINTS; ++i)
+ {
+ B(1, NPOINTS + i) = dNdx(1, i);
+ B(3, i) = dNdx(1, i) / std::sqrt(2);
+ B(3, NPOINTS + i) = dNdx(0, i) / std::sqrt(2);
+ B(0, i) = dNdx(0, i);
}
}
} // detail
@@ -32,47 +33,50 @@ void fillBMatrix2DCartesianPart(DNDX_Type const& dNdx, BMatrixType& b_matrix)
/// Fills a B-matrix based on given shape function dN/dx values.
template <int DisplacementDim,
int NPOINTS,
+ typename BMatrixType,
typename N_Type,
- typename DNDX_Type,
- typename BMatrixType>
-void computeBMatrix(DNDX_Type const& dNdx,
- BMatrixType& b_matrix,
- const bool is_axially_symmetric,
- N_Type const& N,
- const double radius)
+ typename DNDX_Type>
+BMatrixType computeBMatrix(DNDX_Type const& dNdx,
+ N_Type const& N,
+ const double radius,
+ const bool is_axially_symmetric)
{
static_assert(0 < DisplacementDim && DisplacementDim <= 3,
"LinearBMatrix::computeBMatrix: DisplacementDim must be in "
"range [1,3].");
- b_matrix.setZero();
+ BMatrixType B =
+ BMatrixType::Zero(KelvinVectorDimensions<DisplacementDim>::value,
+ NPOINTS * DisplacementDim);
switch (DisplacementDim)
{
case 3:
for (int i = 0; i < NPOINTS; ++i)
{
- b_matrix(2, 2 * NPOINTS + i) = dNdx(2, i);
- b_matrix(4, NPOINTS + i) = dNdx(2, i) / std::sqrt(2);
- b_matrix(4, 2 * NPOINTS + i) = dNdx(1, i) / std::sqrt(2);
- b_matrix(5, i) = dNdx(2, i) / std::sqrt(2);
- b_matrix(5, 2 * NPOINTS + i) = dNdx(0, i) / std::sqrt(2);
+ B(2, 2 * NPOINTS + i) = dNdx(2, i);
+ B(4, NPOINTS + i) = dNdx(2, i) / std::sqrt(2);
+ B(4, 2 * NPOINTS + i) = dNdx(1, i) / std::sqrt(2);
+ B(5, i) = dNdx(2, i) / std::sqrt(2);
+ B(5, 2 * NPOINTS + i) = dNdx(0, i) / std::sqrt(2);
}
- detail::fillBMatrix2DCartesianPart<NPOINTS>(dNdx, b_matrix);
+ detail::fillBMatrix2DCartesianPart<NPOINTS>(dNdx, B);
break;
case 2:
- detail::fillBMatrix2DCartesianPart<NPOINTS>(dNdx, b_matrix);
+ detail::fillBMatrix2DCartesianPart<NPOINTS>(dNdx, B);
if (is_axially_symmetric)
{
for (int i = 0; i < NPOINTS; ++i)
{
- b_matrix(2, i) = N[i] / radius;
+ B(2, i) = N[i] / radius;
}
}
break;
default:
break;
}
+
+ return B;
}
} // namespace LinearBMatrix
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
index 1868619a5d489e5f9f872d0b4a018f0436bf5598..a265823988b4c3020dc98b45197c3daf564ea1f5 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
@@ -45,8 +45,7 @@ struct IntegrationPointData final
// The default generated move-ctor is correctly generated for other
// compilers.
explicit IntegrationPointData(IntegrationPointData&& other)
- : b_matrices(std::move(other.b_matrices)),
- sigma_eff(std::move(other.sigma_eff)),
+ : sigma_eff(std::move(other.sigma_eff)),
sigma_eff_prev(std::move(other.sigma_eff_prev)),
eps(std::move(other.eps)),
eps_prev(std::move(other.eps_prev)),
@@ -59,12 +58,13 @@ struct IntegrationPointData final
typename ShapeMatrixTypeDisplacement::template MatrixType<
DisplacementDim, NPoints * DisplacementDim>
- N_u;
- // typename ShapeMatrixTypeDisplacement::NodalRowVectorType N_u;
- typename BMatricesType::BMatrixType b_matrices;
+ N_u_op;
typename BMatricesType::KelvinVectorType sigma_eff, sigma_eff_prev;
typename BMatricesType::KelvinVectorType eps, eps_prev;
+ typename ShapeMatrixTypeDisplacement::NodalRowVectorType N_u;
+ typename ShapeMatrixTypeDisplacement::GlobalDimNodalMatrixType dNdx_u;
+
typename ShapeMatricesTypePressure::NodalRowVectorType N_p;
typename ShapeMatricesTypePressure::GlobalDimNodalMatrixType dNdx_p;
@@ -88,8 +88,6 @@ struct IntegrationPointData final
double const dt,
DisplacementVectorType const& u)
{
- eps.noalias() = b_matrices * u;
-
auto&& solution = solid_material.integrateStress(
t, x_position, dt, eps_prev, eps, sigma_eff_prev,
*material_state_variables);
@@ -183,12 +181,13 @@ public:
HydroMechanicsLocalAssembler(
MeshLib::Element const& e,
std::size_t const /*local_matrix_size*/,
- bool is_axially_symmetric,
+ bool const is_axially_symmetric,
unsigned const integration_order,
HydroMechanicsProcessData<DisplacementDim>& process_data)
: _process_data(process_data),
_integration_method(integration_order),
- _element(e)
+ _element(e),
+ _is_axially_symmetric(is_axially_symmetric)
{
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
@@ -212,36 +211,27 @@ public:
// displacement (subscript u)
_ip_data.emplace_back(*_process_data.material);
auto& ip_data = _ip_data[ip];
- auto const& sm = shape_matrices_u[ip];
+ auto const& sm_u = shape_matrices_u[ip];
_ip_data[ip].integration_weight =
_integration_method.getWeightedPoint(ip).getWeight() *
- sm.integralMeasure * shape_matrices_u[ip].detJ;
- ip_data.b_matrices.resize(
- kelvin_vector_size,
- ShapeFunctionDisplacement::NPOINTS * DisplacementDim);
-
- auto const x_coord =
- interpolateXCoordinate<ShapeFunctionDisplacement,
- ShapeMatricesTypeDisplacement>(
- e, shape_matrices_u[ip].N);
- LinearBMatrix::computeBMatrix<DisplacementDim,
- ShapeFunctionDisplacement::NPOINTS>(
- shape_matrices_u[ip].dNdx, ip_data.b_matrices,
- is_axially_symmetric, shape_matrices_u[ip].N, x_coord);
+ sm_u.integralMeasure * sm_u.detJ;
ip_data.sigma_eff.resize(kelvin_vector_size);
ip_data.sigma_eff_prev.resize(kelvin_vector_size);
ip_data.eps.resize(kelvin_vector_size);
ip_data.eps_prev.resize(kelvin_vector_size);
- ip_data.N_u = ShapeMatricesTypeDisplacement::template MatrixType<
+ ip_data.N_u_op = ShapeMatricesTypeDisplacement::template MatrixType<
DisplacementDim, displacement_size>::Zero(DisplacementDim,
displacement_size);
for (int i = 0; i < DisplacementDim; ++i)
- ip_data.N_u
+ ip_data.N_u_op
.template block<1, displacement_size / DisplacementDim>(
i, i * displacement_size / DisplacementDim)
- .noalias() = shape_matrices_u[ip].N;
+ .noalias() = sm_u.N;
+
+ ip_data.N_u = sm_u.N;
+ ip_data.dNdx_u = sm_u.dNdx;
ip_data.N_p = shape_matrices_p[ip].N;
ip_data.dNdx_p = shape_matrices_p[ip].dNdx;
@@ -324,11 +314,24 @@ public:
x_position.setIntegrationPoint(ip);
auto const& w = _ip_data[ip].integration_weight;
- auto const& N_p = _ip_data[ip].N_p;
+ auto const& N_u_op = _ip_data[ip].N_u_op;
+
auto const& N_u = _ip_data[ip].N_u;
+ auto const& dNdx_u = _ip_data[ip].dNdx_u;
+
+ auto const& N_p = _ip_data[ip].N_p;
auto const& dNdx_p = _ip_data[ip].dNdx_p;
- auto const& B = _ip_data[ip].b_matrices;
+ auto const x_coord =
+ interpolateXCoordinate<ShapeFunctionDisplacement,
+ ShapeMatricesTypeDisplacement>(_element,
+ N_u);
+ auto const B = LinearBMatrix::computeBMatrix<
+ DisplacementDim, ShapeFunctionDisplacement::NPOINTS,
+ typename BMatricesType::BMatrixType>(dNdx_u, N_u, x_coord,
+ _is_axially_symmetric);
+
+ auto& eps = _ip_data[ip].eps;
auto const& sigma_eff = _ip_data[ip].sigma_eff;
double const S = _process_data.specific_storage(t, x_position)[0];
@@ -341,11 +344,13 @@ public:
auto const porosity = _process_data.porosity(t, x_position)[0];
auto const& b = _process_data.specific_body_force;
auto const& identity2 = MaterialLib::SolidModels::Invariants<
- kelvin_vector_size>::identity2;
+ KelvinVectorDimensions<DisplacementDim>::value>::identity2;
//
// displacement equation, displacement part
//
+ eps.noalias() = B * u;
+
auto C =
_ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u);
@@ -357,7 +362,7 @@ public:
double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
local_rhs.template segment<displacement_size>(displacement_index)
.noalias() -=
- (B.transpose() * sigma_eff - N_u.transpose() * rho * b) * w;
+ (B.transpose() * sigma_eff - N_u_op.transpose() * rho * b) * w;
//
// displacement equation, pressure part
@@ -607,6 +612,7 @@ private:
IntegrationMethod _integration_method;
MeshLib::Element const& _element;
+ bool const _is_axially_symmetric;
SecondaryData<
typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType>
_secondary_data;
diff --git a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h
index d016533b04cc23e3a32b94e6ec34469d64d6208e..b7a6c4c60ad02f774a957dd919ddcff425eea137 100644
--- a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h
+++ b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h
@@ -37,8 +37,9 @@ HydroMechanicsLocalAssemblerFracture<ShapeFunctionDisplacement,
unsigned const integration_order,
HydroMechanicsProcessData<GlobalDim>& process_data)
: HydroMechanicsLocalAssemblerInterface(
- e,
- ShapeFunctionDisplacement::NPOINTS * GlobalDim + ShapeFunctionPressure::NPOINTS,
+ e, is_axially_symmetric,
+ ShapeFunctionDisplacement::NPOINTS * GlobalDim +
+ ShapeFunctionPressure::NPOINTS,
dofIndex_to_localIndex),
_process_data(process_data)
{
diff --git a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerInterface.h b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerInterface.h
index d7dc816721750cb9b21ce1d4be0c0adb1c66f8b0..1524108b4c6c1ea73dfa12c9969f911d44a0934f 100644
--- a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerInterface.h
+++ b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerInterface.h
@@ -33,10 +33,12 @@ class HydroMechanicsLocalAssemblerInterface
{
public:
HydroMechanicsLocalAssemblerInterface(MeshLib::Element const& element,
+ bool const is_axially_symmetric,
std::size_t n_local_size,
std::vector<unsigned>
dofIndex_to_localIndex)
: _element(element),
+ _is_axially_symmetric(is_axially_symmetric),
_dofIndex_to_localIndex(std::move(dofIndex_to_localIndex))
{
_local_u.resize(n_local_size);
@@ -112,6 +114,7 @@ protected:
double const t, Eigen::VectorXd const& local_u) = 0;
MeshLib::Element const& _element;
+ bool const _is_axially_symmetric;
private:
Eigen::VectorXd _local_u;
diff --git a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h
index 9921426409e8d994c93a6d219e54cd4c06c537b2..e2cdc3eca41e577fa220cafc993eea637ccd55b7 100644
--- a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h
+++ b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h
@@ -36,7 +36,7 @@ HydroMechanicsLocalAssemblerMatrix<ShapeFunctionDisplacement,
unsigned const integration_order,
HydroMechanicsProcessData<GlobalDim>& process_data)
: HydroMechanicsLocalAssemblerInterface(
- e,
+ e, is_axially_symmetric,
(n_variables - 1) * ShapeFunctionDisplacement::NPOINTS * GlobalDim +
ShapeFunctionPressure::NPOINTS,
dofIndex_to_localIndex),
@@ -71,21 +71,14 @@ HydroMechanicsLocalAssemblerMatrix<ShapeFunctionDisplacement,
ip_data.integration_weight =
sm_u.detJ * sm_u.integralMeasure *
integration_method.getWeightedPoint(ip).getWeight();
- ip_data.b_matrices.resize(kelvin_vector_size, displacement_size);
ip_data.N_u = sm_u.N;
+ ip_data.dNdx_u = sm_u.dNdx;
ip_data.H_u.setZero(GlobalDim, displacement_size);
for (unsigned i = 0; i < GlobalDim; ++i)
ip_data.H_u.template block<1, displacement_size / GlobalDim>(
i, i * displacement_size / GlobalDim)
.noalias() = ip_data.N_u;
- auto const x_coord =
- interpolateXCoordinate<ShapeFunctionDisplacement,
- ShapeMatricesTypeDisplacement>(e, sm_u.N);
- LinearBMatrix::computeBMatrix<GlobalDim,
- ShapeFunctionDisplacement::NPOINTS>(
- sm_u.dNdx, ip_data.b_matrices, is_axially_symmetric, sm_u.N,
- x_coord);
ip_data.N_p = sm_p.N;
ip_data.dNdx_p = sm_p.dNdx;
@@ -194,11 +187,22 @@ assembleBlockMatricesWithJacobian(
auto& ip_data = _ip_data[ip];
auto const& ip_w = ip_data.integration_weight;
+ auto const& N_u = ip_data.N_u;
+ auto const& dNdx_u = ip_data.dNdx_u;
auto const& N_p = ip_data.N_p;
auto const& dNdx_p = ip_data.dNdx_p;
auto const& H_u = ip_data.H_u;
- auto const& B = ip_data.b_matrices;
+ auto const x_coord =
+ interpolateXCoordinate<ShapeFunctionDisplacement,
+ ShapeMatricesTypeDisplacement>(_element,
+ N_u);
+ auto const B =
+ LinearBMatrix::computeBMatrix<GlobalDim,
+ ShapeFunctionDisplacement::NPOINTS,
+ typename BMatricesType::BMatrixType>(
+ dNdx_u, N_u, x_coord, _is_axially_symmetric);
+
auto const& eps_prev = ip_data.eps_prev;
auto const& sigma_eff_prev = ip_data.sigma_eff_prev;
auto& sigma_eff = ip_data.sigma_eff;
@@ -318,7 +322,6 @@ computeSecondaryVariableConcreteWithBlockVectors(
auto& ip_data = _ip_data[ip];
auto const& dNdx_p = ip_data.dNdx_p;
- auto const& B = ip_data.b_matrices;
auto const& eps_prev = ip_data.eps_prev;
auto const& sigma_eff_prev = ip_data.sigma_eff_prev;
@@ -332,6 +335,19 @@ computeSecondaryVariableConcreteWithBlockVectors(
auto q = ip_data.darcy_velocity.head(GlobalDim);
+ auto const& N_u = ip_data.N_u;
+ auto const& dNdx_u = ip_data.dNdx_u;
+
+ auto const x_coord =
+ interpolateXCoordinate<ShapeFunctionDisplacement,
+ ShapeMatricesTypeDisplacement>(_element,
+ N_u);
+ auto const B =
+ LinearBMatrix::computeBMatrix<GlobalDim,
+ ShapeFunctionDisplacement::NPOINTS,
+ typename BMatricesType::BMatrixType>(
+ dNdx_u, N_u, x_coord, _is_axially_symmetric);
+
eps.noalias() = B * u;
auto&& solution = _ip_data[ip].solid_material.integrateStress(
diff --git a/ProcessLib/LIE/HydroMechanics/LocalAssembler/IntegrationPointDataMatrix.h b/ProcessLib/LIE/HydroMechanics/LocalAssembler/IntegrationPointDataMatrix.h
index dc873c00b9a5688b4bd49cca8dbbf3a0536accb7..e938893427cccae7b6d774df42981010627bc204 100644
--- a/ProcessLib/LIE/HydroMechanics/LocalAssembler/IntegrationPointDataMatrix.h
+++ b/ProcessLib/LIE/HydroMechanics/LocalAssembler/IntegrationPointDataMatrix.h
@@ -36,8 +36,8 @@ struct IntegrationPointDataMatrix final
// compilers.
explicit IntegrationPointDataMatrix(IntegrationPointDataMatrix&& other)
: N_u(std::move(other.N_u)),
+ dNdx_u(std::move(other.dNdx_u)),
H_u(std::move(other.H_u)),
- b_matrices(std::move(other.b_matrices)),
sigma_eff(std::move(other.sigma_eff)),
sigma_eff_prev(std::move(other.sigma_eff_prev)),
eps(std::move(other.eps)),
@@ -53,10 +53,10 @@ struct IntegrationPointDataMatrix final
}
typename ShapeMatrixTypeDisplacement::NodalRowVectorType N_u;
+ typename ShapeMatrixTypeDisplacement::GlobalDimNodalMatrixType dNdx_u;
typename ShapeMatrixTypeDisplacement::template MatrixType<
GlobalDim, NPoints * GlobalDim>
H_u;
- typename BMatricesType::BMatrixType b_matrices;
typename BMatricesType::KelvinVectorType sigma_eff, sigma_eff_prev;
typename BMatricesType::KelvinVectorType eps, eps_prev;
diff --git a/ProcessLib/LIE/SmallDeformation/LocalAssembler/IntegrationPointDataMatrix.h b/ProcessLib/LIE/SmallDeformation/LocalAssembler/IntegrationPointDataMatrix.h
index 73b28dec804b3a13cfda093df13503071e03b367..09760332ac2328c9753355482a969b3abfdd2db3 100644
--- a/ProcessLib/LIE/SmallDeformation/LocalAssembler/IntegrationPointDataMatrix.h
+++ b/ProcessLib/LIE/SmallDeformation/LocalAssembler/IntegrationPointDataMatrix.h
@@ -20,8 +20,8 @@ namespace LIE
{
namespace SmallDeformation
{
-
-template <typename BMatricesType, int DisplacementDim>
+template <typename ShapeMatricesType, typename BMatricesType,
+ int DisplacementDim>
struct IntegrationPointDataMatrix final
{
explicit IntegrationPointDataMatrix(
@@ -36,8 +36,7 @@ struct IntegrationPointDataMatrix final
// The default generated move-ctor is correctly generated for other
// compilers.
explicit IntegrationPointDataMatrix(IntegrationPointDataMatrix&& other)
- : _b_matrices(std::move(other._b_matrices)),
- _sigma(std::move(other._sigma)),
+ : _sigma(std::move(other._sigma)),
_sigma_prev(std::move(other._sigma_prev)),
_eps(std::move(other._eps)),
_eps_prev(std::move(other._eps_prev)),
@@ -49,7 +48,6 @@ struct IntegrationPointDataMatrix final
}
#endif // _MSC_VER
- typename BMatricesType::BMatrixType _b_matrices;
typename BMatricesType::KelvinVectorType _sigma, _sigma_prev;
typename BMatricesType::KelvinVectorType _eps, _eps_prev;
@@ -62,6 +60,9 @@ struct IntegrationPointDataMatrix final
double _detJ;
double _integralMeasure;
+ typename ShapeMatricesType::NodalRowVectorType N;
+ typename ShapeMatricesType::GlobalDimNodalMatrixType dNdx;
+
void pushBackState()
{
_eps_prev = _eps;
diff --git a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix-impl.h b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix-impl.h
index 3418b593a9969cb3c903dd5a82b0cbbb1a8fad59..3c9c657b3e20585a4ae094b00f53ee81d0e6b652 100644
--- a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix-impl.h
+++ b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix-impl.h
@@ -35,20 +35,20 @@ namespace LIE
{
namespace SmallDeformation
{
-
template <typename ShapeFunction, typename IntegrationMethod,
int DisplacementDim>
SmallDeformationLocalAssemblerMatrix<ShapeFunction, IntegrationMethod,
- DisplacementDim>::
+ DisplacementDim>::
SmallDeformationLocalAssemblerMatrix(
MeshLib::Element const& e,
std::size_t const /*local_matrix_size*/,
- bool is_axially_symmetric,
+ bool const is_axially_symmetric,
unsigned const integration_order,
SmallDeformationProcessData<DisplacementDim>& process_data)
: _process_data(process_data),
_integration_method(integration_order),
- _element(e)
+ _element(e),
+ _is_axially_symmetric(is_axially_symmetric)
{
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
@@ -66,20 +66,10 @@ SmallDeformationLocalAssemblerMatrix<ShapeFunction, IntegrationMethod,
_ip_data.emplace_back(*_process_data._material);
auto& ip_data = _ip_data[ip];
auto const& sm = shape_matrices[ip];
+ ip_data.N = sm.N;
+ ip_data.dNdx = sm.dNdx;
ip_data._detJ = sm.detJ;
ip_data._integralMeasure = sm.integralMeasure;
- ip_data._b_matrices.resize(
- KelvinVectorDimensions<DisplacementDim>::value,
- ShapeFunction::NPOINTS * DisplacementDim);
-
- auto const x_coord =
- interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(e,
- sm.N);
- LinearBMatrix::computeBMatrix<DisplacementDim,
- ShapeFunction::NPOINTS>(
- sm.dNdx, ip_data._b_matrices, is_axially_symmetric, sm.N,
- x_coord);
-
ip_data._sigma.resize(KelvinVectorDimensions<DisplacementDim>::value);
ip_data._sigma_prev.resize(
KelvinVectorDimensions<DisplacementDim>::value);
@@ -131,7 +121,17 @@ assembleWithJacobian(
auto const& detJ = _ip_data[ip]._detJ;
auto const& integralMeasure = _ip_data[ip]._integralMeasure;
- auto const& B = _ip_data[ip]._b_matrices;
+ auto const& N = _ip_data[ip].N;
+ auto const& dNdx = _ip_data[ip].dNdx;
+ auto const x_coord =
+ interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(_element,
+ N);
+ auto const B =
+ LinearBMatrix::computeBMatrix<DisplacementDim,
+ ShapeFunction::NPOINTS,
+ typename BMatricesType::BMatrixType>(
+ dNdx, N, x_coord, _is_axially_symmetric);
+
auto const& eps_prev = _ip_data[ip]._eps_prev;
auto const& sigma_prev = _ip_data[ip]._sigma_prev;
diff --git a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix.h b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix.h
index 74217790ea1fb45c066543133c0d74c327ea764e..2bc41ea1e302fec9a8566ceb95c47cd6fa93f42c 100644
--- a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix.h
+++ b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix.h
@@ -42,7 +42,6 @@ public:
using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
using BMatricesType = BMatrixPolicyType<ShapeFunction, DisplacementDim>;
- using BMatrixType = typename BMatricesType::BMatrixType;
using StiffnessMatrixType = typename BMatricesType::StiffnessMatrixType;
using NodalForceVectorType = typename BMatricesType::NodalForceVectorType;
using NodalDisplacementVectorType =
@@ -159,13 +158,15 @@ private:
SmallDeformationProcessData<DisplacementDim>& _process_data;
- std::vector<IntegrationPointDataMatrix<BMatricesType, DisplacementDim>,
- Eigen::aligned_allocator<
- IntegrationPointDataMatrix<BMatricesType, DisplacementDim>>>
+ std::vector<IntegrationPointDataMatrix<ShapeMatricesType, BMatricesType,
+ DisplacementDim>,
+ Eigen::aligned_allocator<IntegrationPointDataMatrix<
+ ShapeMatricesType, BMatricesType, DisplacementDim>>>
_ip_data;
IntegrationMethod _integration_method;
MeshLib::Element const& _element;
+ bool const _is_axially_symmetric;
SecondaryData<typename ShapeMatrices::ShapeType> _secondary_data;
};
diff --git a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture-impl.h b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture-impl.h
index 84b70ec46453964f98a6e7ab64a2e7cbe6f852a1..3ff08e674ecc3a95714cacebc57f9d35399ef5ec 100644
--- a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture-impl.h
+++ b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture-impl.h
@@ -42,11 +42,12 @@ namespace LIE
{
namespace SmallDeformation
{
-
-template <typename ShapeFunction, typename IntegrationMethod,
+template <typename ShapeFunction,
+ typename IntegrationMethod,
int DisplacementDim>
-SmallDeformationLocalAssemblerMatrixNearFracture<ShapeFunction, IntegrationMethod,
- DisplacementDim>::
+SmallDeformationLocalAssemblerMatrixNearFracture<ShapeFunction,
+ IntegrationMethod,
+ DisplacementDim>::
SmallDeformationLocalAssemblerMatrixNearFracture(
MeshLib::Element const& e,
std::size_t const n_variables,
@@ -60,12 +61,18 @@ SmallDeformationLocalAssemblerMatrixNearFracture<ShapeFunction, IntegrationMetho
dofIndex_to_localIndex),
_process_data(process_data),
_integration_method(integration_order),
- _shape_matrices(
- initShapeMatrices<ShapeFunction, ShapeMatricesType,
- IntegrationMethod, DisplacementDim>(
- e, is_axially_symmetric, _integration_method)),
- _element(e)
+ _element(e),
+ _is_axially_symmetric(is_axially_symmetric)
{
+ std::vector<
+ ShapeMatrices,
+ Eigen::aligned_allocator<typename ShapeMatricesType::ShapeMatrices>>
+ shape_matrices = initShapeMatrices<ShapeFunction,
+ ShapeMatricesType,
+ IntegrationMethod,
+ DisplacementDim>(
+ e, is_axially_symmetric, _integration_method);
+
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
@@ -76,20 +83,11 @@ SmallDeformationLocalAssemblerMatrixNearFracture<ShapeFunction, IntegrationMetho
{
_ip_data.emplace_back(*_process_data._material);
auto& ip_data = _ip_data[ip];
- auto const& sm = _shape_matrices[ip];
+ auto const& sm = shape_matrices[ip];
+ ip_data.N = sm.N;
+ ip_data.dNdx = sm.dNdx;
ip_data._detJ = sm.detJ;
ip_data._integralMeasure = sm.integralMeasure;
- ip_data._b_matrices.resize(
- KelvinVectorDimensions<DisplacementDim>::value,
- ShapeFunction::NPOINTS * DisplacementDim);
-
- auto const x_coord =
- interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(e,
- sm.N);
- LinearBMatrix::computeBMatrix<DisplacementDim,
- ShapeFunction::NPOINTS>(
- sm.dNdx, ip_data._b_matrices, is_axially_symmetric, sm.N,
- x_coord);
ip_data._sigma.resize(KelvinVectorDimensions<DisplacementDim>::value);
ip_data._sigma_prev.resize(
@@ -184,13 +182,15 @@ assembleWithJacobian(
{
x_position.setIntegrationPoint(ip);
- auto const& sm = _shape_matrices[ip];
auto &ip_data = _ip_data[ip];
auto const& wp = _integration_method.getWeightedPoint(ip);
auto const ip_factor = ip_data._detJ * wp.getWeight() * ip_data._integralMeasure;
+ auto const& N = ip_data.N;
+ auto const& dNdx = ip_data.dNdx;
+
// levelset functions
- auto const ip_physical_coords = computePhysicalCoordinates(_element, sm.N);
+ auto const ip_physical_coords = computePhysicalCoordinates(_element, N);
std::vector<double> levelsets(n_fractures);
for (unsigned i = 0; i < n_fractures; i++)
levelsets[i] = calculateLevelSetFunction(*_fracture_props[i],
@@ -201,8 +201,16 @@ assembleWithJacobian(
for (unsigned i=0; i<n_fractures; i++)
nodal_total_u += levelsets[i] * vec_nodal_g[i];
+ auto const x_coord =
+ interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(_element,
+ N);
+ auto const B =
+ LinearBMatrix::computeBMatrix<DisplacementDim,
+ ShapeFunction::NPOINTS,
+ typename BMatricesType::BMatrixType>(
+ dNdx, N, x_coord, _is_axially_symmetric);
+
// strain, stress
- auto const& B = ip_data._b_matrices;
auto const& eps_prev = ip_data._eps_prev;
auto const& sigma_prev = ip_data._sigma_prev;
diff --git a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture.h b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture.h
index 82ba949ef6805bd39671beff3fa15aff986404d7..095bfbb89b3e5a00861f7c16f2eb58c2829b0ea4 100644
--- a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture.h
+++ b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture.h
@@ -161,17 +161,16 @@ private:
SmallDeformationProcessData<DisplacementDim>& _process_data;
std::vector<FractureProperty*> _fracture_props;
- std::vector<IntegrationPointDataMatrix<BMatricesType, DisplacementDim>,
- Eigen::aligned_allocator<
- IntegrationPointDataMatrix<BMatricesType, DisplacementDim>>>
+ std::vector<IntegrationPointDataMatrix<ShapeMatricesType, BMatricesType,
+ DisplacementDim>,
+ Eigen::aligned_allocator<IntegrationPointDataMatrix<
+ ShapeMatricesType, BMatricesType, DisplacementDim>>>
_ip_data;
IntegrationMethod _integration_method;
- std::vector<ShapeMatrices, Eigen::aligned_allocator<
- typename ShapeMatricesType::ShapeMatrices>>
- _shape_matrices;
MeshLib::Element const& _element;
+ bool const _is_axially_symmetric;
SecondaryData<typename ShapeMatrices::ShapeType> _secondary_data;
};
diff --git a/ProcessLib/SmallDeformation/SmallDeformationFEM.h b/ProcessLib/SmallDeformation/SmallDeformationFEM.h
index e8493f15fce106d4ba992576c66eccd72ae319d1..d0e1d6d1a55f71c747285b11f6ad06e4a8be3022 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationFEM.h
+++ b/ProcessLib/SmallDeformation/SmallDeformationFEM.h
@@ -31,7 +31,8 @@ namespace ProcessLib
{
namespace SmallDeformation
{
-template <typename BMatricesType, int DisplacementDim>
+template <typename BMatricesType, typename ShapeMatricesType,
+ int DisplacementDim>
struct IntegrationPointData final
{
explicit IntegrationPointData(
@@ -46,7 +47,6 @@ struct IntegrationPointData final
// The default generated move-ctor is correctly generated for other
// compilers.
explicit IntegrationPointData(IntegrationPointData&& other)
- : b_matrices(std::move(other.b_matrices)),
sigma(std::move(other.sigma)),
sigma_prev(std::move(other.sigma_prev)),
eps(std::move(other.eps)),
@@ -58,7 +58,6 @@ struct IntegrationPointData final
}
#endif // _MSC_VER
- typename BMatricesType::BMatrixType b_matrices;
typename BMatricesType::KelvinVectorType sigma, sigma_prev;
typename BMatricesType::KelvinVectorType eps, eps_prev;
@@ -68,6 +67,8 @@ struct IntegrationPointData final
material_state_variables;
double integration_weight;
+ typename ShapeMatricesType::NodalRowVectorType N;
+ typename ShapeMatricesType::GlobalDimNodalMatrixType dNdx;
void pushBackState()
{
@@ -153,12 +154,13 @@ public:
SmallDeformationLocalAssembler(
MeshLib::Element const& e,
std::size_t const /*local_matrix_size*/,
- bool is_axially_symmetric,
+ bool const is_axially_symmetric,
unsigned const integration_order,
SmallDeformationProcessData<DisplacementDim>& process_data)
: _process_data(process_data),
_integration_method(integration_order),
- _element(e)
+ _element(e),
+ _is_axially_symmetric(is_axially_symmetric)
{
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
@@ -179,25 +181,15 @@ public:
_ip_data[ip].integration_weight =
_integration_method.getWeightedPoint(ip).getWeight() *
sm.integralMeasure * sm.detJ;
- ip_data.b_matrices.resize(
- KelvinVectorDimensions<DisplacementDim>::value,
- ShapeFunction::NPOINTS * DisplacementDim);
- auto const x_coord =
- interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(e,
- sm.N);
- LinearBMatrix::computeBMatrix<DisplacementDim,
- ShapeFunction::NPOINTS>(
- sm.dNdx, ip_data.b_matrices, is_axially_symmetric, sm.N,
- x_coord);
-
- ip_data.sigma.resize(
- KelvinVectorDimensions<DisplacementDim>::value);
+ ip_data.N = sm.N;
+ ip_data.dNdx = sm.dNdx;
+
+ ip_data.sigma.resize(KelvinVectorDimensions<DisplacementDim>::value);
ip_data.sigma_prev.resize(
KelvinVectorDimensions<DisplacementDim>::value);
ip_data.eps.resize(KelvinVectorDimensions<DisplacementDim>::value);
- ip_data.eps_prev.resize(
- KelvinVectorDimensions<DisplacementDim>::value);
+ ip_data.eps_prev.resize(KelvinVectorDimensions<DisplacementDim>::value);
_secondary_data.N[ip] = shape_matrices[ip].N;
}
@@ -240,8 +232,17 @@ public:
{
x_position.setIntegrationPoint(ip);
auto const& w = _ip_data[ip].integration_weight;
+ auto const& N = _ip_data[ip].N;
+ auto const& dNdx = _ip_data[ip].dNdx;
+
+ auto const x_coord =
+ interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
+ _element, N);
+ auto const B = LinearBMatrix::computeBMatrix<
+ DisplacementDim, ShapeFunction::NPOINTS,
+ typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
+ _is_axially_symmetric);
- auto const& B = _ip_data[ip].b_matrices;
auto const& eps_prev = _ip_data[ip].eps_prev;
auto const& sigma_prev = _ip_data[ip].sigma_prev;
@@ -286,9 +287,10 @@ public:
std::vector<double>& nodal_values) const override
{
return ProcessLib::SmallDeformation::getNodalForces<
- DisplacementDim, ShapeFunction::NPOINTS,
- NodalDisplacementVectorType>(nodal_values, _integration_method,
- _ip_data, _element.getID());
+ DisplacementDim, ShapeFunction, ShapeMatricesType,
+ NodalDisplacementVectorType, typename BMatricesType::BMatrixType>(
+ nodal_values, _integration_method, _ip_data, _element,
+ _is_axially_symmetric);
}
Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
@@ -413,14 +415,16 @@ private:
SmallDeformationProcessData<DisplacementDim>& _process_data;
- std::vector<IntegrationPointData<BMatricesType, DisplacementDim>,
- Eigen::aligned_allocator<
- IntegrationPointData<BMatricesType, DisplacementDim>>>
+ std::vector<
+ IntegrationPointData<BMatricesType, ShapeMatricesType, DisplacementDim>,
+ Eigen::aligned_allocator<IntegrationPointData<
+ BMatricesType, ShapeMatricesType, DisplacementDim>>>
_ip_data;
IntegrationMethod _integration_method;
MeshLib::Element const& _element;
SecondaryData<typename ShapeMatrices::ShapeType> _secondary_data;
+ bool const _is_axially_symmetric;
};
template <typename ShapeFunction, typename IntegrationMethod,
diff --git a/ProcessLib/SmallDeformationCommon/Common.h b/ProcessLib/SmallDeformationCommon/Common.h
index 7787f381dcb20598c0582c8ed96a4ee1b39a7647..5b65b7e31aed95cda56b1238d5d61fb1d11ba1f4 100644
--- a/ProcessLib/SmallDeformationCommon/Common.h
+++ b/ProcessLib/SmallDeformationCommon/Common.h
@@ -16,6 +16,9 @@
#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
#include "NumLib/DOF/DOFTableUtil.h"
+#include "ProcessLib/Deformation/BMatrixPolicy.h"
+#include "ProcessLib/Deformation/LinearBMatrix.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
namespace ProcessLib
{
@@ -29,30 +32,37 @@ struct NodalForceCalculationInterface
virtual ~NodalForceCalculationInterface() = default;
};
-template <int DisplacementDim, int NPoints,
- typename NodalDisplacementVectorType, typename IPData,
- typename IntegrationMethod>
+template <int DisplacementDim, typename ShapeFunction,
+ typename ShapeMatricesType, typename NodalDisplacementVectorType,
+ typename BMatrixType, typename IPData, typename IntegrationMethod>
std::vector<double> const& getNodalForces(
std::vector<double>& nodal_values,
IntegrationMethod const& _integration_method, IPData const& _ip_data,
- int element_id)
+ MeshLib::Element const& element, bool const is_axially_symmetric)
{
nodal_values.clear();
auto local_b = MathLib::createZeroedVector<NodalDisplacementVectorType>(
- nodal_values, NPoints * DisplacementDim);
+ nodal_values, ShapeFunction::NPOINTS * DisplacementDim);
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
SpatialPosition x_position;
- x_position.setElementID(element_id);
+ x_position.setElementID(element.getID());
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
x_position.setIntegrationPoint(ip);
auto const& w = _ip_data[ip].integration_weight;
- auto const& B = _ip_data[ip].b_matrices;
+ auto const x_coord =
+ interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
+ element, _ip_data[ip].N);
+ auto const B =
+ LinearBMatrix::computeBMatrix<DisplacementDim,
+ ShapeFunction::NPOINTS, BMatrixType>(
+ _ip_data[ip].dNdx, _ip_data[ip].N, x_coord,
+ is_axially_symmetric);
auto& sigma = _ip_data[ip].sigma;
local_b.noalias() += B.transpose() * sigma * w;