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