swh:1:snp:f521c49ab17ef7db6ec70b2430e1ed203f50383f
Tip revision: 10a1660b266d5ce990a32ce7098825b2d898b899 authored by montoyav on 22 March 2021, 16:59:16 UTC
Change of nomenclature (in the documentation) of the physical constants in the Material Library
Change of nomenclature (in the documentation) of the physical constants in the Material Library
Tip revision: 10a1660
NormalTractionBoundaryConditionLocalAssembler.h
/**
* \file
* \copyright
* Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*
*/
#pragma once
#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
#include "MeshLib/Elements/FaceRule.h"
#include "NumLib/DOF/DOFTableUtil.h"
#include "ParameterLib/Parameter.h"
#include "GenericNaturalBoundaryConditionLocalAssembler.h"
namespace ProcessLib
{
namespace NormalTractionBoundaryCondition
{
template <typename ShapeMatricesType>
struct IntegrationPointData final
{
IntegrationPointData(
typename ShapeMatricesType::ShapeMatrices::ShapeType const N_,
typename ShapeMatricesType::GlobalDimVectorType const n_,
double const integration_weight_)
: N(N_), n(n_), integration_weight(integration_weight_)
{
}
typename ShapeMatricesType::ShapeMatrices::ShapeType const N;
typename ShapeMatricesType::GlobalDimVectorType const n;
double const integration_weight;
};
class NormalTractionBoundaryConditionLocalAssemblerInterface
{
public:
virtual void assemble(
std::size_t const id,
NumLib::LocalToGlobalIndexMap const& dof_table_boundary, double const t,
std::vector<GlobalVector*> const& /*x*/, GlobalMatrix& /*K*/,
GlobalVector& b, GlobalMatrix* /*Jac*/) = 0;
virtual ~NormalTractionBoundaryConditionLocalAssemblerInterface() = default;
};
template <typename ShapeFunctionDisplacement, typename IntegrationMethod,
unsigned GlobalDim>
class NormalTractionBoundaryConditionLocalAssembler final
: public NormalTractionBoundaryConditionLocalAssemblerInterface
{
public:
using ShapeMatricesType =
ShapeMatrixPolicyType<ShapeFunctionDisplacement, GlobalDim>;
using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType;
using NodalVectorType = typename ShapeMatricesType::NodalVectorType;
NormalTractionBoundaryConditionLocalAssembler(
MeshLib::Element const& e,
std::size_t const local_matrix_size,
bool const is_axially_symmetric,
unsigned const integration_order,
ParameterLib::Parameter<double> const& pressure)
: _integration_method(integration_order),
_pressure(pressure),
_element(e)
{
_local_rhs.setZero(local_matrix_size);
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
_ip_data.reserve(n_integration_points);
auto const shape_matrices_u =
NumLib::initShapeMatrices<ShapeFunctionDisplacement,
ShapeMatricesType, GlobalDim>(
e, is_axially_symmetric, _integration_method);
GlobalDimVectorType element_normal(GlobalDim);
// TODO Extend to rotated 2d meshes and line elements.
if (e.getGeomType() == MeshLib::MeshElemType::LINE)
{
Eigen::Vector3d const v1 =
Eigen::Map<Eigen::Vector3d const>(e.getNode(1)->getCoords()) -
Eigen::Map<Eigen::Vector3d const>(e.getNode(0)->getCoords());
element_normal[0] = -v1[1];
element_normal[1] = v1[0];
element_normal.normalize();
}
else
{
auto const n = MeshLib::FaceRule::getSurfaceNormal(&e).normalized();
for (decltype(GlobalDim) i = 0; i < GlobalDim; ++i)
{
element_normal[i] = n[i];
}
}
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
double const integration_weight =
_integration_method.getWeightedPoint(ip).getWeight() *
shape_matrices_u[ip].integralMeasure *
shape_matrices_u[ip].detJ;
_ip_data.emplace_back(shape_matrices_u[ip].N, element_normal,
integration_weight);
}
}
void assemble(std::size_t const id,
NumLib::LocalToGlobalIndexMap const& dof_table_boundary,
double const t, std::vector<GlobalVector*> const& /*x*/,
GlobalMatrix& /*K*/, GlobalVector& local_rhs,
GlobalMatrix* /*Jac*/) override
{
_local_rhs.setZero();
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
NodalVectorType pressure =
_pressure.getNodalValuesOnElement(_element, t);
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
auto const& w = _ip_data[ip].integration_weight;
auto const& N = _ip_data[ip].N;
auto const& n = _ip_data[ip].n;
typename ShapeMatricesType::template MatrixType<GlobalDim,
displacement_size>
N_u = ShapeMatricesType::template MatrixType<
GlobalDim, displacement_size>::Zero(GlobalDim,
displacement_size);
for (int i = 0; i < static_cast<int>(GlobalDim); ++i)
{
N_u.template block<1, displacement_size / GlobalDim>(
i, i * displacement_size / GlobalDim)
.noalias() = N;
}
_local_rhs.noalias() -= n.transpose() * N_u * pressure.dot(N) * w;
}
auto const indices = NumLib::getIndices(id, dof_table_boundary);
local_rhs.add(indices, _local_rhs);
}
private:
IntegrationMethod const _integration_method;
ParameterLib::Parameter<double> const& _pressure;
static const int displacement_size =
ShapeFunctionDisplacement::NPOINTS * GlobalDim;
std::vector<
IntegrationPointData<ShapeMatricesType>,
Eigen::aligned_allocator<IntegrationPointData<ShapeMatricesType>>>
_ip_data;
typename ShapeMatricesType::template VectorType<displacement_size>
_local_rhs;
MeshLib::Element const& _element;
public:
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
};
} // namespace NormalTractionBoundaryCondition
} // namespace ProcessLib