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