OGS
RobinBoundaryConditionLocalAssembler.h
Go to the documentation of this file.
1 
11 #pragma once
12 
15 #include "ParameterLib/Parameter.h"
16 
17 namespace ProcessLib
18 {
20 {
24 };
25 
26 template <typename ShapeFunction, typename IntegrationMethod, int GlobalDim>
29  ShapeFunction, IntegrationMethod, GlobalDim>
30 {
32  ShapeFunction, IntegrationMethod, GlobalDim>;
34 
35 public:
37  std::size_t const local_matrix_size,
38  bool is_axially_symmetric,
39  unsigned const integration_order,
40  RobinBoundaryConditionData const& data)
41  : Base(e, is_axially_symmetric, integration_order),
42  _data(data),
43  _local_K(local_matrix_size, local_matrix_size),
44  _local_rhs(local_matrix_size)
45  {
46  }
47 
48  // TODO also implement derivative for Jacobian in Newton scheme.
49  void assemble(std::size_t const id,
50  NumLib::LocalToGlobalIndexMap const& dof_table_boundary,
51  double const t, std::vector<GlobalVector*> const& /*x*/,
52  int const /*process_id*/, GlobalMatrix& K, GlobalVector& b,
53  GlobalMatrix* /*Jac*/) override
54  {
55  _local_K.setZero();
56  _local_rhs.setZero();
57 
58  unsigned const n_integration_points =
59  Base::_integration_method.getNumberOfPoints();
60 
61  typename Base::NodalVectorType const alpha =
63  .template topRows<ShapeFunction::MeshElement::n_all_nodes>();
64  typename Base::NodalVectorType const u_0 =
66  .template topRows<ShapeFunction::MeshElement::n_all_nodes>();
67 
68  for (unsigned ip = 0; ip < n_integration_points; ++ip)
69  {
70  auto const& ip_data = Base::_ns_and_weights[ip];
71  auto const& N = ip_data.N;
72  auto const& w = ip_data.weight;
73 
74  ParameterLib::SpatialPosition const position{
75  std::nullopt, Base::_element.getID(), ip,
77  NumLib::interpolateCoordinates<ShapeFunction,
79  Base::_element, N))};
80 
81  double integral_measure = 1.0;
83  {
84  integral_measure = (*_data.integral_measure)(t, position)[0];
85  }
86 
87  // flux = alpha * ( u_0 - u )
88  // adding a alpha term to the diagonal of the stiffness matrix
89  // and a alpha * u_0 term to the rhs vector
90  _local_K.diagonal().noalias() +=
91  N * alpha.dot(N) * w * integral_measure;
92  _local_rhs.noalias() +=
93  N * alpha.dot(N) * u_0.dot(N) * w * integral_measure;
94  }
95 
96  auto const indices = NumLib::getIndices(id, dof_table_boundary);
98  _local_K);
99  b.add(indices, _local_rhs);
100  }
101 
102 private:
104 
107 
108 public:
110 };
111 
112 } // namespace ProcessLib
int add(IndexType row, IndexType col, double val)
Definition: EigenMatrix.h:87
Global vector based on Eigen vector.
Definition: EigenVector.h:26
void add(IndexType rowId, double v)
add entry
Definition: EigenVector.h:80
virtual std::size_t getID() const final
Returns the ID of the element.
Definition: Element.h:82
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 &, int const, GlobalMatrix &K, GlobalVector &b, GlobalMatrix *) override
RobinBoundaryConditionLocalAssembler(MeshLib::Element const &e, std::size_t const local_matrix_size, bool is_axially_symmetric, unsigned const integration_order, RobinBoundaryConditionData const &data)
MathLib::TemplatePoint< double, 3 > Point3d
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:163
ParameterLib::Parameter< double > const *const integral_measure
ParameterLib::Parameter< double > const & alpha
ParameterLib::Parameter< double > const & u_0