OGS
RobinBoundaryConditionLocalAssembler.h
Go to the documentation of this file.
1
11#pragma once
12
17
18namespace ProcessLib
19{
26
27template <typename ShapeFunction, int GlobalDim>
30 GlobalDim>
31{
32 using Base =
35
36public:
38 MeshLib::Element const& e,
39 std::size_t const local_matrix_size,
40 NumLib::GenericIntegrationMethod const& integration_method,
41 bool is_axially_symmetric,
43 : Base(e, is_axially_symmetric, integration_method),
44 _data(data),
45 _local_K(local_matrix_size, local_matrix_size),
46 _local_rhs(local_matrix_size)
47 {
48 }
49
50 void assemble(std::size_t const id,
51 NumLib::LocalToGlobalIndexMap const& dof_table_boundary,
52 double const t, std::vector<GlobalVector*> const& xs,
53 int const process_id, GlobalMatrix& K, GlobalVector& b,
54 GlobalMatrix* Jac) override
55 {
56 _local_K.setZero();
57 _local_rhs.setZero();
58
59 auto& x = *xs[process_id];
60
61 auto const indices = NumLib::getIndices(id, dof_table_boundary);
62 auto const local_x = x.get(indices);
63 auto const u =
64 MathLib::toVector<Eigen::Matrix<double, ShapeFunction::NPOINTS, 1>>(
65 local_x, ShapeFunction::NPOINTS);
66
67 unsigned const n_integration_points =
69
70 typename Base::NodalVectorType const alpha =
72 .template topRows<ShapeFunction::MeshElement::n_all_nodes>();
73 typename Base::NodalVectorType const u_0 =
75 .template topRows<ShapeFunction::MeshElement::n_all_nodes>();
76
77 for (unsigned ip = 0; ip < n_integration_points; ++ip)
78 {
79 auto const& ip_data = Base::_ns_and_weights[ip];
80 auto const& N = ip_data.N;
81 auto const& w = ip_data.weight;
82
83 ParameterLib::SpatialPosition const position{
84 std::nullopt, Base::_element.getID(), ip,
88 Base::_element, N))};
89
90 double integral_measure = 1.0;
92 {
93 integral_measure = (*_data.integral_measure)(t, position)[0];
94 }
95
96 double const a = alpha.dot(N) * w * integral_measure;
97
98 // The local K matrix is used for both, the Newton and Picard
99 // methods.
100 _local_K.noalias() += N.transpose() * N * a;
101
102 if (Jac != nullptr)
103 {
104 _local_rhs.noalias() -= N.transpose() * (u - u_0).dot(N) * a;
105 }
106 else
107 {
108 _local_rhs.noalias() += N.transpose() * (u_0.dot(N) * a);
109 }
110 }
111
112 b.add(indices, _local_rhs);
113 if (Jac != nullptr)
114 {
116 indices),
117 _local_K);
118 }
119 else
120 {
122 indices),
123 _local_K);
124 }
125 }
126
127private:
129
132
133public:
135};
136
137} // namespace ProcessLib
int add(IndexType row, IndexType col, double val)
Definition EigenMatrix.h:86
Global vector based on Eigen vector.
Definition EigenVector.h:25
void add(IndexType rowId, double v)
add entry
Definition EigenVector.h:76
std::size_t getID() const
Returns the ID of the element.
Definition Element.h:89
std::vector< NAndWeight, Eigen::aligned_allocator< NAndWeight > > const _ns_and_weights
void assemble(std::size_t const id, NumLib::LocalToGlobalIndexMap const &dof_table_boundary, double const t, std::vector< GlobalVector * > const &xs, int const process_id, GlobalMatrix &K, GlobalVector &b, GlobalMatrix *Jac) override
RobinBoundaryConditionLocalAssembler(MeshLib::Element const &e, std::size_t const local_matrix_size, NumLib::GenericIntegrationMethod const &integration_method, bool is_axially_symmetric, RobinBoundaryConditionData const &data)
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
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
ParameterLib::Parameter< double > const *const integral_measure
ParameterLib::Parameter< double > const & alpha
ParameterLib::Parameter< double > const & u_0