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