OGS
NormalTractionBoundaryConditionLocalAssembler.h
Go to the documentation of this file.
1
11#pragma once
12
18
19namespace ProcessLib
20{
21namespace NormalTractionBoundaryCondition
22{
23template <typename ShapeMatricesType>
25{
27 typename ShapeMatricesType::ShapeMatrices::ShapeType const N_,
28 typename ShapeMatricesType::GlobalDimVectorType const n_,
29 double const integration_weight_)
30 : N(N_), n(n_), integration_weight(integration_weight_)
31 {
32 }
33
34 typename ShapeMatricesType::ShapeMatrices::ShapeType const N;
35 typename ShapeMatricesType::GlobalDimVectorType const n;
36 double const integration_weight;
37};
38
40{
41public:
42 virtual void assemble(
43 std::size_t const id,
44 NumLib::LocalToGlobalIndexMap const& dof_table_boundary, double const t,
45 std::vector<GlobalVector*> const& /*x*/, GlobalMatrix* /*K*/,
46 GlobalVector& b, GlobalMatrix* /*Jac*/) = 0;
48};
49
50template <typename ShapeFunctionDisplacement, int GlobalDim>
53{
54public:
59
61 MeshLib::Element const& e,
62 std::size_t const local_matrix_size,
63 NumLib::GenericIntegrationMethod const& integration_method,
64 bool const is_axially_symmetric,
65 ParameterLib::Parameter<double> const& pressure)
66 : _integration_method(integration_method),
67 _pressure(pressure),
68 _element(e)
69 {
70 _local_rhs.setZero(local_matrix_size);
71
72 unsigned const n_integration_points =
74
75 _ip_data.reserve(n_integration_points);
76
77 auto const shape_matrices_u =
78 NumLib::initShapeMatrices<ShapeFunctionDisplacement,
79 ShapeMatricesType, GlobalDim>(
80 e, is_axially_symmetric, _integration_method);
81
82 GlobalDimVectorType element_normal(GlobalDim);
83
84 // TODO Extend to rotated 2d meshes and line elements.
86 {
87 Eigen::Vector3d const v1 = e.getNode(1)->asEigenVector3d() -
89 element_normal[GlobalDim - 1] = 0;
90 element_normal[0] = -v1[1];
91 element_normal[1] = v1[0];
92 element_normal.normalize();
93 }
94 else
95 {
96 auto const n = MeshLib::FaceRule::getSurfaceNormal(e).normalized();
97 for (int i = 0; i < GlobalDim; ++i)
98 {
99 element_normal[i] = n[i];
100 }
101 }
102
103 for (unsigned ip = 0; ip < n_integration_points; ip++)
104 {
105 double const integration_weight =
107 shape_matrices_u[ip].integralMeasure *
108 shape_matrices_u[ip].detJ;
109
110 _ip_data.emplace_back(shape_matrices_u[ip].N, element_normal,
111 integration_weight);
112 }
113 }
114
115 void assemble(std::size_t const id,
116 NumLib::LocalToGlobalIndexMap const& dof_table_boundary,
117 double const t, std::vector<GlobalVector*> const& /*x*/,
118 GlobalMatrix* /*K*/, GlobalVector& local_rhs,
119 GlobalMatrix* /*Jac*/) override
120 {
121 _local_rhs.setZero();
122
123 unsigned const n_integration_points =
125
126 NodalVectorType pressure =
128
129 for (unsigned ip = 0; ip < n_integration_points; ip++)
130 {
131 auto const& w = _ip_data[ip].integration_weight;
132 auto const& N = _ip_data[ip].N;
133 auto const& n = _ip_data[ip].n;
134
135 typename ShapeMatricesType::template MatrixType<GlobalDim,
137 N_u = ShapeMatricesType::template MatrixType<
138 GlobalDim, displacement_size>::Zero(GlobalDim,
140 for (int i = 0; i < GlobalDim; ++i)
141 {
142 N_u.template block<1, displacement_size / GlobalDim>(
143 i, i * displacement_size / GlobalDim)
144 .noalias() = N;
145 }
146
147 _local_rhs.noalias() -= n.transpose() * N_u * pressure.dot(N) * w;
148 }
149
150 auto const indices = NumLib::getIndices(id, dof_table_boundary);
151 local_rhs.add(indices, _local_rhs);
152 }
153
154private:
157
158 static const int displacement_size =
159 ShapeFunctionDisplacement::NPOINTS * GlobalDim;
160 std::vector<
162 Eigen::aligned_allocator<IntegrationPointData<ShapeMatricesType>>>
164
165 typename ShapeMatricesType::template VectorType<displacement_size>
167
169
170public:
172};
173
174} // namespace NormalTractionBoundaryCondition
175} // namespace ProcessLib
Global vector based on Eigen vector.
Definition EigenVector.h:25
void add(IndexType rowId, double v)
add entry
Definition EigenVector.h:76
Eigen::Vector3d const & asEigenVector3d() const
Definition Point3d.h:64
double getWeight() const
virtual MeshElemType getGeomType() const =0
virtual const Node * getNode(unsigned idx) const =0
static Eigen::Vector3d getSurfaceNormal(Element const &e)
Returns the surface normal of a 2D element.
Definition FaceRule.cpp:40
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
virtual void assemble(std::size_t const id, NumLib::LocalToGlobalIndexMap const &dof_table_boundary, double const t, std::vector< GlobalVector * > const &, GlobalMatrix *, GlobalVector &b, GlobalMatrix *)=0
void assemble(std::size_t const id, NumLib::LocalToGlobalIndexMap const &dof_table_boundary, double const t, std::vector< GlobalVector * > const &, GlobalMatrix *, GlobalVector &local_rhs, GlobalMatrix *) override
NormalTractionBoundaryConditionLocalAssembler(MeshLib::Element const &e, std::size_t const local_matrix_size, NumLib::GenericIntegrationMethod const &integration_method, bool const is_axially_symmetric, ParameterLib::Parameter< double > const &pressure)
std::vector< IntegrationPointData< ShapeMatricesType >, Eigen::aligned_allocator< IntegrationPointData< ShapeMatricesType > > > _ip_data
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, IntegrationMethod const &integration_method)
VectorType< GlobalDim > GlobalDimVectorType
VectorType< ShapeFunction::NPOINTS > NodalVectorType
virtual Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > getNodalValuesOnElement(MeshLib::Element const &element, double const t) const
Returns a matrix of values for all nodes of the given element.
Definition Parameter.h:164
IntegrationPointData(typename ShapeMatricesType::ShapeMatrices::ShapeType const N_, typename ShapeMatricesType::GlobalDimVectorType const n_, double const integration_weight_)