OGS
BcAndStLocalAssemblerImpl.h
Go to the documentation of this file.
1
11#pragma once
12
15#include "NsAndWeight.h"
19
21{
22
31template <typename BcOrStData, typename ShapeFunction,
32 typename LowerOrderShapeFunction, typename IntegrationMethod,
33 int GlobalDim>
35{
36public:
37 using Traits =
39
40private:
42
43public:
45 bool is_axially_symmetric,
46 unsigned const integration_order,
47 BcOrStData const& data)
48 : bc_or_st_data(data),
49 element(e),
50 integration_method(integration_order),
52 computeNsAndWeights<ShapeFunction, LowerOrderShapeFunction,
53 GlobalDim>(e, is_axially_symmetric,
55 {
56 }
57
60 void assemble(std::size_t const boundary_element_id,
61 NumLib::LocalToGlobalIndexMap const& dof_table_boundary,
62 double const t, GlobalVector const& x, GlobalVector& b,
63 GlobalMatrix* const Jac) const
64 {
65 // Variables for input data
66 Eigen::MatrixXd const primary_variables_mat = ProcessLib::
69 dof_table_boundary, x);
70
71 // Variables for output data
72 auto const& indices_this_component =
73 dof_table_boundary(boundary_element_id,
75 .rows;
76 std::vector<GlobalIndexType> indices_all_components_for_Jac;
77
78 auto const num_dof_this_component = indices_this_component.size();
79
80 Eigen::VectorXd local_rhs =
81 Eigen::VectorXd::Zero(num_dof_this_component);
82 Eigen::MatrixXd local_Jac;
83
84 if (Jac)
85 {
86 indices_all_components_for_Jac =
87 NumLib::getIndices(boundary_element_id, dof_table_boundary);
88
89 auto const num_dof_all_components =
90 indices_all_components_for_Jac.size();
91
92 local_Jac = Eigen::MatrixXd::Zero(num_dof_this_component,
93 num_dof_all_components);
94 }
95
96 // Variables for temporary data
97 auto const num_comp_total =
98 dof_table_boundary.getNumberOfGlobalComponents();
99 std::vector<double> prim_vars_data(num_comp_total);
100 auto prim_vars = MathLib::toVector(prim_vars_data);
101
102 // Variables for the integration loop
103 unsigned const num_integration_points =
104 integration_method.getNumberOfPoints();
106 ShapeFunction, ShapeMatrixPolicy>(element);
107
108 for (unsigned ip = 0; ip < num_integration_points; ip++)
109 {
110 auto const& ns_and_weight = nss_and_weights[ip];
111 auto const coords = interpolateCoords(ns_and_weight, fe);
112
114 primary_variables_mat,
116 nss_and_weights[ip], prim_vars);
117
118 auto const [flag, flux, dFlux] =
119 bc_or_st_data.getFlagAndFluxAndDFlux(t, coords, prim_vars_data);
120
121 if (!flag)
122 {
123 // No flux value for this integration point. Skip assembly of
124 // the entire element.
125 return;
126 }
127
128 assembleLocalRhs(flux, ns_and_weight, local_rhs);
129
130 if (static_cast<int>(dFlux.size()) != num_comp_total)
131 {
132 // This strict check is technically mandatory only if a Jacobian
133 // is assembled. However, it is done as a consistency check also
134 // for cases without Jacobian assembly.
135 OGS_FATAL(
136 "The Python BC or ST must return the derivative of the "
137 "flux w.r.t. each component of each primary variable. "
138 "{:d} components expected. {:d} components returned from "
139 "Python."
140 "The expected number of components depends on multiple "
141 "factors. In the case of vectorial primary variables (such "
142 "as displacement), it will depend on the space dimension. "
143 "In the case of staggered coupling, the number of "
144 "components can be different for each staggered process "
145 "and different from monolithic coupling.",
146 num_comp_total, dFlux.size());
147 }
148
149 if (Jac)
150 {
151 assembleLocalJacobian(dFlux, ns_and_weight, local_Jac);
152 }
153 }
154
155 b.add(indices_this_component, local_rhs);
156
157 if (Jac)
158 {
160 indices_this_component, indices_all_components_for_Jac};
161 Jac->add(rci, local_Jac);
162 }
163 }
164
165private:
171 std::array<double, 3> interpolateCoords(
172 typename Traits::NsAndWeight const& ns_and_weight,
174 fe) const
175 {
176 auto const& N_higher = ns_and_weight.NHigherOrder();
177 return fe.interpolateCoordinates(N_higher);
178 }
179
182 void assembleLocalRhs(double const flux,
183 typename Traits::NsAndWeight const& ns_and_weight,
184 Eigen::VectorXd& local_rhs) const
185 {
186 auto const w = ns_and_weight.weight();
187 auto const N = ns_and_weight.N(bc_or_st_data.shape_function_order);
188
189 local_rhs.noalias() += N.transpose() * (flux * w);
190 }
191
205 std::vector<double> const& dFlux,
206 typename Traits::NsAndWeight const& ns_and_weight,
207 Eigen::MatrixXd& local_Jac) const
208 {
209 auto const& pv_refs =
211 auto const weight = ns_and_weight.weight();
212
213 auto const this_shp_fct_order = bc_or_st_data.shape_function_order;
214 auto const N_this = ns_and_weight.N(this_shp_fct_order);
215
216 Eigen::Index left = 0;
217
218 for (auto pv_ref : pv_refs)
219 {
220 auto const& pv = pv_ref.get();
221 auto const other_num_comp = pv.getNumberOfGlobalComponents();
222 auto const other_shp_fct_order = pv.getShapeFunctionOrder();
223 auto const N_other = ns_and_weight.N(other_shp_fct_order);
224 auto const width = N_other.size();
225
226 for (decltype(+other_num_comp) other_comp = 0;
227 other_comp < other_num_comp;
228 ++other_comp)
229 {
230 // The assignment -= takes into account the sign convention
231 // of 1st-order in time in ODE systems in OpenGeoSys.
232 local_Jac.middleCols(left, width).noalias() -=
233 N_this.transpose() * (dFlux[other_comp] * weight) * N_other;
234
235 left += width;
236 }
237 }
238
239 assert(left == local_Jac.cols());
240 }
241
242public:
245
246 IntegrationMethod const integration_method;
247 std::vector<typename Traits::NsAndWeight> const nss_and_weights;
248};
249
250} // namespace ProcessLib::BoundaryConditionAndSourceTerm::Python
#define OGS_FATAL(...)
Definition: Error.h:26
int add(IndexType row, IndexType col, double val)
Definition: EigenMatrix.h:87
Global vector based on Eigen vector.
Definition: EigenVector.h:28
void add(IndexType rowId, double v)
add entry
Definition: EigenVector.h:79
std::size_t getID() const
Get id of the mesh.
Definition: Mesh.h:115
Template class for isoparametric elements.
std::array< double, 3 > interpolateCoordinates(typename ShapeMatrices::ShapeType const &N) const
Interpolates the coordinates of the element with the given shape matrix.
Eigen::Map< const Vector > toVector(std::vector< double > const &data, Eigen::VectorXd::Index size)
Creates an Eigen mapped vector from the given data vector.
static const double t
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
NumLib::TemplateIsoparametric< ShapeFunction, ShapeMatricesType > createIsoparametricFiniteElement(MeshLib::Element const &e)
auto computeNsAndWeights(MeshLib::Element const &element, bool const is_axially_symmetric, IntegrationMethod const &integration_method)
Definition: NsAndWeight.h:143
Eigen::MatrixXd collectDofsToMatrix(MeshLib::Element const &element, std::size_t const mesh_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x)
void interpolate(Eigen::MatrixXd const &primary_variables_mat, std::vector< std::reference_wrapper< ProcessVariable > > const &pv_refs, NsAndWeight const &ns_and_weight, Eigen::Ref< Eigen::VectorXd > interpolated_primary_variables)
void assembleLocalJacobian(std::vector< double > const &dFlux, typename Traits::NsAndWeight const &ns_and_weight, Eigen::MatrixXd &local_Jac) const
std::array< double, 3 > interpolateCoords(typename Traits::NsAndWeight const &ns_and_weight, NumLib::TemplateIsoparametric< ShapeFunction, ShapeMatrixPolicy > const &fe) const
BcAndStLocalAssemblerImpl(MeshLib::Element const &e, bool is_axially_symmetric, unsigned const integration_order, BcOrStData const &data)
void assembleLocalRhs(double const flux, typename Traits::NsAndWeight const &ns_and_weight, Eigen::VectorXd &local_rhs) const
void assemble(std::size_t const boundary_element_id, NumLib::LocalToGlobalIndexMap const &dof_table_boundary, double const t, GlobalVector const &x, GlobalVector &b, GlobalMatrix *const Jac) const
MeshLib::Mesh const & bc_or_st_mesh
The domain, where this BC or ST will be applied.
Definition: BcOrStData.h:46
std::vector< std::reference_wrapper< ProcessVariable > > const & all_process_variables_for_this_process
Definition: BcOrStData.h:51
Eigen::Ref< const Eigen::RowVectorXd > N(unsigned order) const
Definition: NsAndWeight.h:51
Collects common type aliases needed when working with NsAndWeight.
Definition: NsAndWeight.h:119
ShapeMatrixPolicyType< ShapeFunction, GlobalDim > ShapeMatrixPolicy
Definition: NsAndWeight.h:120