OGS 6.2.0-97-g4a610c866
PythonBoundaryConditionLocalAssembler.h
Go to the documentation of this file.
1 
10 #pragma once
11 
13 
17 
19 
20 namespace ProcessLib
21 {
25 {
26 };
27 
28 template <typename ShapeFunction, typename IntegrationMethod,
29  unsigned GlobalDim>
32  ShapeFunction, IntegrationMethod, GlobalDim>
33 {
35  ShapeFunction, IntegrationMethod, GlobalDim>;
36 
37 public:
39  MeshLib::Element const& e,
40  std::size_t const /*local_matrix_size*/,
41  bool is_axially_symmetric,
42  unsigned const integration_order,
43  PythonBoundaryConditionData const& data)
44  : Base(e, is_axially_symmetric, integration_order), _data(data)
45  {
46  }
47 
48  void assemble(std::size_t const boundary_element_id,
49  NumLib::LocalToGlobalIndexMap const& dof_table_boundary,
50  double const t, const GlobalVector& x, GlobalMatrix& /*K*/,
51  GlobalVector& b, GlobalMatrix* Jac) override
52  {
53  using ShapeMatricesType =
55  using FemType =
57  FemType fe(*static_cast<const typename ShapeFunction::MeshElement*>(
58  &this->_element));
59 
60  unsigned const num_integration_points =
61  Base::_integration_method.getNumberOfPoints();
62  auto const num_var = _data.dof_table_bulk.getNumberOfVariables();
63  auto const num_nodes = Base::_element.getNumberOfNodes();
64  auto const num_comp_total =
65  _data.dof_table_bulk.getNumberOfComponents();
66 
67  auto const& bulk_node_ids_map =
68  *_data.boundary_mesh.getProperties()
69  .template getPropertyVector<std::size_t>(
70  "bulk_node_ids", MeshLib::MeshItemType::Node, 1);
71 
72  // gather primary variables
73  Eigen::MatrixXd primary_variables_mat(num_nodes, num_comp_total);
74  for (int var = 0; var < num_var; ++var)
75  {
76  auto const num_comp =
77  _data.dof_table_bulk.getNumberOfVariableComponents(var);
78  for (int comp = 0; comp < num_comp; ++comp)
79  {
80  auto const global_component =
81  dof_table_boundary.getGlobalComponent(var, comp);
82 
83  for (unsigned element_node_id = 0; element_node_id < num_nodes;
84  ++element_node_id)
85  {
86  auto const* const node =
87  Base::_element.getNode(element_node_id);
88  auto const boundary_node_id = node->getID();
89  auto const bulk_node_id =
90  bulk_node_ids_map[boundary_node_id];
91  MeshLib::Location loc{_data.bulk_mesh_id,
93  bulk_node_id};
94  auto const dof_idx =
95  _data.dof_table_bulk.getGlobalIndex(loc, var, comp);
96  if (dof_idx == NumLib::MeshComponentMap::nop)
97  {
98  // TODO extend Python BC to mixed FEM ansatz functions
99  OGS_FATAL(
100  "No d.o.f. found for (node=%d, var=%d, comp=%d). "
101  "That might be due to the use of mixed FEM ansatz "
102  "functions, which is currently not supported by "
103  "the implementation of Python BCs. That excludes, "
104  "e.g., the HM process.",
105  bulk_node_id, var, comp);
106  }
107  primary_variables_mat(element_node_id, global_component) =
108  x[dof_idx];
109  }
110  }
111  }
112 
113  Eigen::VectorXd local_rhs = Eigen::VectorXd::Zero(num_nodes);
114  Eigen::MatrixXd local_Jac =
115  Eigen::MatrixXd::Zero(num_nodes, num_nodes * num_comp_total);
116  std::vector<double> prim_vars_data(num_comp_total);
117  auto prim_vars = MathLib::toVector(prim_vars_data);
118 
119  for (unsigned ip = 0; ip < num_integration_points; ip++)
120  {
121  auto const& N = Base::_ns_and_weights[ip].N;
122  auto const& w = Base::_ns_and_weights[ip].weight;
123  auto const coords = fe.interpolateCoordinates(N);
124  prim_vars =
125  N * primary_variables_mat; // Assumption: all primary variables
126  // have same shape functions.
127  auto const flag_flux_dFlux =
128  _data.bc_object->getFlux(t, coords, prim_vars_data);
129  if (!_data.bc_object->isOverriddenNatural())
130  {
131  // getFlux() is not overridden in Python, so we can skip the
132  // whole BC assembly (i.e., for all boundary elements).
134  }
135 
136  if (!std::get<0>(flag_flux_dFlux))
137  {
138  // No flux value for this integration point. Skip assembly of
139  // the entire element.
140  return;
141  }
142  auto const flux = std::get<1>(flag_flux_dFlux);
143  auto const& dFlux = std::get<2>(flag_flux_dFlux);
144 
145  local_rhs.noalias() += N * (flux * w);
146 
147  if (static_cast<int>(dFlux.size()) != num_comp_total)
148  {
149  // This strict check is technically mandatory only if a Jacobian
150  // is assembled. However, it is done as a consistency check also
151  // for cases without Jacobian assembly.
152  OGS_FATAL(
153  "The Python BC must return the derivative of the flux "
154  "w.r.t. each primary variable. %d components expected. %d "
155  "components returned from Python.",
156  num_comp_total, dFlux.size());
157  }
158 
159  if (Jac)
160  {
161  for (int comp = 0; comp < num_comp_total; ++comp)
162  {
163  auto const top = 0;
164  auto const left = comp * num_nodes;
165  auto const width = num_nodes;
166  auto const height = num_nodes;
167  // The assignement -= takes into account the sign convention
168  // of 1st-order in time ODE systems in OpenGeoSys.
169  local_Jac.block(top, left, width, height).noalias() -=
170  N.transpose() * (dFlux[comp] * w) * N;
171  }
172  }
173  }
174 
175  auto const& indices_specific_component =
176  dof_table_boundary(boundary_element_id, _data.global_component_id)
177  .rows;
178  b.add(indices_specific_component, local_rhs);
179 
180  if (Jac)
181  {
182  // only assemble a block of the Jacobian, not the whole local matrix
183  auto const indices_all_components =
184  NumLib::getIndices(boundary_element_id, dof_table_boundary);
186  indices_specific_component, indices_all_components};
187  Jac->add(rci, local_Jac);
188  }
189  }
190 
191 private:
193 };
194 
195 } // namespace ProcessLib
int getGlobalComponent(int const variable_id, int const component_id) const
void assemble(std::size_t const boundary_element_id, NumLib::LocalToGlobalIndexMap const &dof_table_boundary, double const t, const GlobalVector &x, GlobalMatrix &, GlobalVector &b, GlobalMatrix *Jac) override
Eigen::Map< const Vector > toVector(std::vector< double > const &data, Eigen::VectorXd::Index size)
Creates an Eigen mapped vector from the given data vector.
PythonBoundaryConditionLocalAssembler(MeshLib::Element const &e, std::size_t const, bool is_axially_symmetric, unsigned const integration_order, PythonBoundaryConditionData const &data)
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
#define OGS_FATAL(fmt,...)
Definition: Error.h:63
Template class for isoparametric elements.
static NUMLIB_EXPORT GlobalIndexType const nop