OGS 6.2.1-97-g73d1aeda3
PythonBoundaryCondition.cpp
Go to the documentation of this file.
1 
11 
12 #include <pybind11/pybind11.h>
13 #include <iostream>
14 
18 
19 namespace
20 {
25 {
26 public:
28  explicit FlushStdoutGuard(bool const flush) : _flush(flush)
29  {
30  if (!flush)
31  {
32  return;
33  }
34 
35  LOGOG_COUT << std::flush;
36  }
37 
40  {
41  if (!_flush)
42  {
43  return;
44  }
45 
46  using namespace pybind11::literals;
47  pybind11::print("end"_a = "", "flush"_a = true);
48  }
49 
50 private:
52  const bool _flush;
53 };
54 } // anonymous namespace
55 
56 namespace ProcessLib
57 {
58 PythonBoundaryCondition::PythonBoundaryCondition(
60  unsigned const integration_order,
61  unsigned const shapefunction_order,
62  unsigned const global_dim,
63  bool const flush_stdout)
64  : _bc_data(std::move(bc_data)), _flush_stdout(flush_stdout)
65 {
66  std::vector<MeshLib::Node*> const& bc_nodes =
68  MeshLib::MeshSubset bc_mesh_subset(_bc_data.boundary_mesh, bc_nodes);
69 
70  // Create local DOF table from the bc mesh subset for the given variable and
71  // component id.
73  std::move(bc_mesh_subset));
74 
75  createLocalAssemblers<PythonBoundaryConditionLocalAssembler>(
77  shapefunction_order, _local_assemblers,
78  _bc_data.boundary_mesh.isAxiallySymmetric(), integration_order,
79  _bc_data);
80 }
81 
83  const double t, GlobalVector const& x,
85 {
86  FlushStdoutGuard guard(_flush_stdout);
87  (void)guard;
88 
89  auto const nodes = _bc_data.boundary_mesh.getNodes();
90 
91  auto const& bulk_node_ids_map =
93  "bulk_node_ids", MeshLib::MeshItemType::Node, 1);
94 
95  bc_values.ids.clear();
96  bc_values.values.clear();
97 
98  bc_values.ids.reserve(_bc_data.boundary_mesh.getNumberOfNodes());
99  bc_values.values.reserve(_bc_data.boundary_mesh.getNumberOfNodes());
100 
101  std::vector<double> primary_variables;
102 
103  for (auto const* node : _bc_data.boundary_mesh.getNodes())
104  {
105  auto const boundary_node_id = node->getID();
106  auto const bulk_node_id = bulk_node_ids_map[boundary_node_id];
107 
108  // gather primary variable values
109  primary_variables.clear();
110  auto const num_var = _dof_table_boundary->getNumberOfVariables();
111  for (int var = 0; var < num_var; ++var)
112  {
113  auto const num_comp =
114  _dof_table_boundary->getNumberOfVariableComponents(var);
115  for (int comp = 0; comp < num_comp; ++comp)
116  {
119  bulk_node_id};
120  auto const dof_idx =
121  _bc_data.dof_table_bulk.getGlobalIndex(loc, var, comp);
122 
123  if (dof_idx == NumLib::MeshComponentMap::nop)
124  {
125  // TODO extend Python BC to mixed FEM ansatz functions
126  OGS_FATAL(
127  "No d.o.f. found for (node=%d, var=%d, comp=%d). "
128  "That might be due to the use of mixed FEM ansatz "
129  "functions, which is currently not supported by "
130  "the implementation of Python BCs. That excludes, "
131  "e.g., the HM process.",
132  bulk_node_id, var, comp);
133  }
134 
135  primary_variables.push_back(x[dof_idx]);
136  }
137  }
138 
139  auto* xs = nodes[boundary_node_id]->getCoords(); // TODO DDC problems?
140  auto pair_flag_value = _bc_data.bc_object->getDirichletBCValue(
141  t, {xs[0], xs[1], xs[2]}, boundary_node_id, primary_variables);
143  {
144  DBUG(
145  "Method `getDirichletBCValue' not overridden in Python "
146  "script.");
147  return;
148  }
149 
150  if (!pair_flag_value.first)
151  {
152  continue;
153  }
154 
155  MeshLib::Location l(_bc_data.bulk_mesh_id, MeshLib::MeshItemType::Node,
156  bulk_node_id);
157  const auto dof_idx = _bc_data.dof_table_bulk.getGlobalIndex(
159  if (dof_idx == NumLib::MeshComponentMap::nop)
160  {
161  OGS_FATAL(
162  "Logic error. This error should already have occured while "
163  "gathering primary variables. Something nasty is going on!");
164  }
165 
166  // For the DDC approach (e.g. with PETSc option), the negative
167  // index of g_idx means that the entry by that index is a ghost
168  // one, which should be dropped. Especially for PETSc routines
169  // MatZeroRows and MatZeroRowsColumns, which are called to apply
170  // the Dirichlet BC, the negative index is not accepted like
171  // other matrix or vector PETSc routines. Therefore, the
172  // following if-condition is applied.
173  if (dof_idx >= 0)
174  {
175  bc_values.ids.emplace_back(dof_idx);
176  bc_values.values.emplace_back(pair_flag_value.second);
177  }
178  }
179 }
180 
182  const GlobalVector& x,
183  GlobalMatrix& K, GlobalVector& b,
184  GlobalMatrix* Jac)
185 {
186  FlushStdoutGuard guard(_flush_stdout);
187 
188  try
189  {
192  _local_assemblers, *_dof_table_boundary, t, x, K, b, Jac);
193  }
194  catch (MethodNotOverriddenInDerivedClassException const& /*e*/)
195  {
196  DBUG("Method `getFlux' not overridden in Python script.");
197  }
198 }
199 
200 std::unique_ptr<PythonBoundaryCondition> createPythonBoundaryCondition(
201  BaseLib::ConfigTree const& config, MeshLib::Mesh const& boundary_mesh,
202  NumLib::LocalToGlobalIndexMap const& dof_table, std::size_t bulk_mesh_id,
203  int const variable_id, int const component_id,
204  unsigned const integration_order, unsigned const shapefunction_order,
205  unsigned const global_dim)
206 {
208  config.checkConfigParameter("type", "Python");
209 
211  auto const bc_object = config.getConfigParameter<std::string>("bc_object");
213  auto const flush_stdout = config.getConfigParameter("flush_stdout", false);
214 
215  // Evaluate Python code in scope of main module
216  pybind11::object scope =
217  pybind11::module::import("__main__").attr("__dict__");
218 
219  if (!scope.contains(bc_object))
220  {
221  OGS_FATAL(
222  "Function `%s' is not defined in the python script file, or there "
223  "was no python script file specified.",
224  bc_object.c_str());
225  }
226 
227  auto* bc = scope[bc_object.c_str()]
229 
230  if (variable_id >= static_cast<int>(dof_table.getNumberOfVariables()) ||
231  component_id >= dof_table.getNumberOfVariableComponents(variable_id))
232  {
233  OGS_FATAL(
234  "Variable id or component id too high. Actual values: (%d, %d), "
235  "maximum values: (%d, %d).",
236  variable_id, component_id, dof_table.getNumberOfVariables(),
237  dof_table.getNumberOfVariableComponents(variable_id));
238  }
239 
240  // In case of partitioned mesh the boundary could be empty, i.e. there is no
241  // boundary condition.
242 #ifdef USE_PETSC
243  // This can be extracted to createBoundaryCondition() but then the config
244  // parameters are not read and will cause an error.
245  // TODO (naumov): Add a function to ConfigTree for skipping the tags of the
246  // subtree and move the code up in createBoundaryCondition().
247  if (boundary_mesh.getDimension() == 0 &&
248  boundary_mesh.getNumberOfNodes() == 0 &&
249  boundary_mesh.getNumberOfElements() == 0)
250  {
251  return nullptr;
252  }
253 #endif // USE_PETSC
254 
255  return std::make_unique<PythonBoundaryCondition>(
257  bc, dof_table, bulk_mesh_id,
258  dof_table.getGlobalComponent(variable_id, component_id),
259  boundary_mesh},
260  integration_order, shapefunction_order, global_dim, flush_stdout);
261 }
262 
263 } // namespace ProcessLib
int getGlobalComponent(int const variable_id, int const component_id) const
std::size_t getNumberOfNodes() const
Get the number of nodes.
Definition: Mesh.h:99
PythonBoundaryConditionPythonSideInterface * bc_object
Python object computing BC values.
FlushStdoutGuard(bool const flush)
Optionally flushes C++ stdout before running Python code.
NumLib::LocalToGlobalIndexMap const & dof_table_bulk
DOF table of the entire domain.
std::unique_ptr< PythonBoundaryCondition > createPythonBoundaryCondition(BaseLib::ConfigTree const &config, MeshLib::Mesh const &boundary_mesh, NumLib::LocalToGlobalIndexMap const &dof_table, std::size_t bulk_mesh_id, int const variable_id, int const component_id, unsigned const integration_order, unsigned const shapefunction_order, unsigned const global_dim)
Creates a new PythonBoundaryCondition object.
PythonBoundaryConditionData _bc_data
Auxiliary data.
std::size_t const bulk_mesh_id
Mesh ID of the entire domain.
static void executeMemberOnDereferenced(Method method, Container const &container, Args &&... args)
int getNumberOfVariableComponents(int variable_id) const
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
Definition: Mesh.h:105
std::vector< std::unique_ptr< GenericNaturalBoundaryConditionLocalAssemblerInterface > > _local_assemblers
Local assemblers for all elements of the boundary mesh.
MeshLib::Properties & getProperties()
Definition: Mesh.h:134
bool isAxiallySymmetric() const
Definition: Mesh.h:137
T getConfigParameter(std::string const &param) const
std::unique_ptr< NumLib::LocalToGlobalIndexMap > _dof_table_boundary
Local dof table for the boundary mesh.
~FlushStdoutGuard()
Optionally flushes Python&#39;s stdout after running Python code.
void applyNaturalBC(const double t, const GlobalVector &x, GlobalMatrix &K, GlobalVector &b, GlobalMatrix *Jac) override
virtual void assemble(std::size_t const id, NumLib::LocalToGlobalIndexMap const &dof_table_boundary, double const t, const GlobalVector &x, GlobalMatrix &K, GlobalVector &b, GlobalMatrix *Jac)=0
const MeshLib::Mesh & boundary_mesh
The boundary mesh, i.e., the domain of this BC.
void checkConfigParameter(std::string const &param, T const &value) const
void getEssentialBCValues(const double t, const GlobalVector &x, NumLib::IndexValueVector< GlobalIndexType > &bc_values) const override
Writes the values of essential BCs to bc_values.
std::vector< IndexType > ids
GlobalIndexType getGlobalIndex(MeshLib::Location const &l, int const variable_id, int const component_id) const
PropertyVector< T > const * getPropertyVector(std::string const &name) const
Definition: Properties.h:119
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition: Mesh.h:108
virtual std::pair< bool, double > getDirichletBCValue(double, std::array< double, 3 >, std::size_t, std::vector< double > const &) const
LocalToGlobalIndexMap * deriveBoundaryConstrainedMap(int const variable_id, std::vector< int > const &component_ids, MeshLib::MeshSubset &&new_mesh_subset) const
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
Definition: Mesh.h:81
std::size_t getNumberOfElements() const
Get the number of elements.
Definition: Mesh.h:96
#define OGS_FATAL(fmt,...)
Definition: Error.h:63
A subset of nodes on a single mesh.
Definition: MeshSubset.h:26
std::vector< double > values
static NUMLIB_EXPORT GlobalIndexType const nop