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