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