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