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"
23
24namespace
25{
27 std::size_t const nnodes)
28{
29 bc_values.ids.clear();
30 bc_values.values.clear();
31
32 bc_values.ids.reserve(nnodes);
33 bc_values.values.reserve(nnodes);
34}
35
37 NumLib::LocalToGlobalIndexMap const& dof_table,
38 std::vector<std::reference_wrapper<ProcessLib::ProcessVariable>> const& pvs)
39{
40 auto const num_vars_dt = dof_table.getNumberOfVariables();
41 auto const num_vars_pv = pvs.size();
42
43 if (static_cast<std::size_t>(num_vars_dt) != num_vars_pv)
44 {
46 "The number of variables in the d.o.f. table does not match the "
47 "number of process variables: {} != {}.",
48 num_vars_dt, num_vars_pv);
49 }
50
51 for (std::size_t var = 0; var < num_vars_pv; ++var)
52 {
53 auto const num_comp_dt = dof_table.getNumberOfVariableComponents(var);
54 auto const num_comp_pv = pvs[var].get().getNumberOfGlobalComponents();
55
56 if (num_comp_dt != num_comp_pv)
57 {
59 "The number of components of variable #{} in the d.o.f. table "
60 "does not match the number of components of process variable "
61 "#{} ({}): {} != {}.",
62 var, var, pvs[var].get().getName(), num_vars_dt, num_vars_pv);
63 }
64 }
65}
66} // anonymous namespace
67
68namespace ProcessLib
69{
71 PythonBcData&& bc_data, unsigned const integration_order,
72 bool const flush_stdout, unsigned const bulk_mesh_dimension,
73 NumLib::LocalToGlobalIndexMap const& dof_table_bulk)
74 : _bc_data(std::move(bc_data)), _flush_stdout(flush_stdout)
75{
76 checkConsistency(dof_table_bulk,
78
79 std::vector<MeshLib::Node*> const& bc_nodes =
81 MeshLib::MeshSubset bc_mesh_subset(_bc_data.bc_or_st_mesh, bc_nodes);
82
84 dof_table_bulk.deriveBoundaryConstrainedMap(std::move(bc_mesh_subset));
85
88 bulk_mesh_dimension, _bc_data.bc_or_st_mesh.getElements(),
90 NumLib::IntegrationOrder{integration_order},
92}
93
95 const double t, GlobalVector const& x,
97{
98 FlushStdoutGuard guard(_flush_stdout);
99 (void)guard;
100
101 auto const& boundary_mesh = _bc_data.bc_or_st_mesh;
102 auto const& boundary_nodes = boundary_mesh.getNodes();
103 auto const* bc_object = _bc_data.bc_or_st_object;
104
105 initBCValues(bc_values, boundary_nodes.size());
106
107 std::vector<double> primary_variables;
108
109 for (auto const* boundary_node : boundary_nodes)
110 {
111 auto const boundary_node_id = boundary_node->getID();
112 auto const dof_idx = getDofIdx(boundary_node_id);
113
114 if (dof_idx == NumLib::MeshComponentMap::nop)
115 {
116 // This d.o.f. has not been found. This can be the case, e.g., for
117 // Taylor-Hood elements, where the lower order field has no d.o.f.
118 // on higher order nodes.
119 continue;
120 }
121
122 if (dof_idx < 0)
123 {
124 // For the DDC approach (e.g. with PETSc option) a negative
125 // index means that this entry is a ghost entry and should be
126 // dropped.
127 continue;
128 }
129
130 collectPrimaryVariables(primary_variables, *boundary_node, x);
131
132 auto const [apply_bc, bc_value] = bc_object->getDirichletBCValue(
133 t, {(*boundary_node)[0], (*boundary_node)[1], (*boundary_node)[2]},
134 boundary_node_id, primary_variables);
135
136 if (!bc_object->isOverriddenEssential())
137 {
138 DBUG(
139 "Method `getDirichletBCValue' not overridden in Python "
140 "script.");
141 return;
142 }
143
144 if (!apply_bc)
145 {
146 continue;
147 }
148
149 bc_values.ids.emplace_back(dof_idx);
150 bc_values.values.emplace_back(bc_value);
151 }
152}
153
155 std::size_t const boundary_node_id) const
156{
158 MeshLib::MeshItemType::Node, boundary_node_id};
159 return _dof_table_boundary->getGlobalIndex(loc,
161}
162
164 std::size_t const boundary_node_id, int const var, int const comp) const
165{
167 MeshLib::MeshItemType::Node, boundary_node_id};
168 return _dof_table_boundary->getGlobalIndex(loc, var, comp);
169}
170
172 std::vector<double>& primary_variables, MeshLib::Node const& boundary_node,
173 GlobalVector const& x) const
174{
175 primary_variables.clear();
176 auto const num_var = _dof_table_boundary->getNumberOfVariables();
177 auto const boundary_node_id = boundary_node.getID();
178
179 for (int var = 0; var < num_var; ++var)
180 {
181 auto const num_comp =
182 _dof_table_boundary->getNumberOfVariableComponents(var);
183 for (int comp = 0; comp < num_comp; ++comp)
184 {
185 auto const dof_idx = getDofIdx(boundary_node_id, var, comp);
186
187 double const pv_value =
189 ? x[dof_idx]
190 : interpolateToHigherOrderNode(x, var, comp, boundary_node);
191
192 primary_variables.push_back(pv_value);
193 }
194 }
195}
196
198 GlobalVector const& x, int const var, int const comp,
199 MeshLib::Node const& boundary_node) const
200{
201 auto const& boundary_elements =
203 boundary_node.getID());
204
205 if (boundary_elements.size() != 1)
206 {
207 DBUG(
208 "Boundary node {} is associated with {} elements in the boundary "
209 "mesh.",
210 boundary_node.getID(),
211 boundary_elements.size());
212 }
213
214 // Interpolations on all associated elements should return the same. Just
215 // pick any of them.
216 auto const& boundary_element = *boundary_elements.front();
217
218 assert(boundary_element.getNumberOfBaseNodes() <
219 boundary_element.getNumberOfNodes() &&
220 "We expect that the boundary element is a higher order element. "
221 "Otherwise no interpolation should take place.");
222
223 // Search local node id
224 auto const local_node_id_within_boundary_element =
225 getNodeIDinElement(boundary_element, &boundary_node);
226
227 auto const boundary_element_id = boundary_element.getID();
228
229 // Assumption: all boundary elements have a local assembler (in the same
230 // order)
231 auto const& loc_asm = *_local_assemblers[boundary_element_id];
232
233 return loc_asm.interpolate(local_node_id_within_boundary_element,
234 *_dof_table_boundary, x, var, comp);
235}
236
238 const double t, std::vector<GlobalVector*> const& x, int const process_id,
240{
241 FlushStdoutGuard guard(_flush_stdout);
242
243 try
244 {
247 _local_assemblers, *_dof_table_boundary, t, x, process_id, K, b,
248 Jac);
249 }
251 {
252 DBUG("Method `getFlux' not overridden in Python script.");
253 }
254}
255
256std::unique_ptr<PythonBoundaryCondition> createPythonBoundaryCondition(
257 BaseLib::ConfigTree const& config, MeshLib::Mesh const& boundary_mesh,
258 NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
259 MeshLib::Mesh const& bulk_mesh, int const variable_id,
260 int const component_id, unsigned const integration_order,
261 unsigned const shapefunction_order,
262 std::vector<std::reference_wrapper<ProcessVariable>> const&
263 all_process_variables_for_this_process)
264{
266 config.checkConfigParameter("type", "Python");
267
269 auto const bc_object = config.getConfigParameter<std::string>("bc_object");
271 auto const flush_stdout = config.getConfigParameter("flush_stdout", false);
272
273 // Evaluate Python code in scope of main module
274 pybind11::object scope =
275 pybind11::module::import("__main__").attr("__dict__");
276
277 if (!scope.contains(bc_object))
278 {
279 OGS_FATAL(
280 "Function `{:s}' is not defined in the python script file, or "
281 "there was no python script file specified.",
282 bc_object);
283 }
284
285 auto* bc = scope[bc_object.c_str()]
287
288 if (variable_id >=
289 static_cast<int>(dof_table_bulk.getNumberOfVariables()) ||
290 component_id >=
291 dof_table_bulk.getNumberOfVariableComponents(variable_id))
292 {
293 OGS_FATAL(
294 "Variable id or component id too high. Actual values: ({:d}, "
295 "{:d}), maximum values: ({:d}, {:d}).",
296 variable_id, component_id, dof_table_bulk.getNumberOfVariables(),
297 dof_table_bulk.getNumberOfVariableComponents(variable_id));
298 }
299
300 // In case of partitioned mesh the boundary could be empty, i.e. there is no
301 // boundary condition.
302#ifdef USE_PETSC
303 // This can be extracted to createBoundaryCondition() but then the config
304 // parameters are not read and will cause an error.
305 // TODO (naumov): Add a function to ConfigTree for skipping the tags of the
306 // subtree and move the code up in createBoundaryCondition().
307 if (boundary_mesh.getDimension() == 0 &&
308 boundary_mesh.getNumberOfNodes() == 0 &&
309 boundary_mesh.getNumberOfElements() == 0)
310 {
311 return nullptr;
312 }
313#endif // USE_PETSC
314
315 return std::make_unique<PythonBoundaryCondition>(
317 {bc, dof_table_bulk.getGlobalComponent(variable_id, component_id),
318 boundary_mesh, all_process_variables_for_this_process,
319 shapefunction_order}},
320 integration_order, flush_stdout, bulk_mesh.getDimension(),
321 dof_table_bulk);
322}
323
324} // namespace ProcessLib
#define OGS_FATAL(...)
Definition: Error.h:26
GlobalMatrix::IndexType GlobalIndexType
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition: Logging.h:30
std::string getName(std::string const &line)
Returns the name/title from the "Zone"-description.
T getConfigParameter(std::string const &param) const
void checkConfigParameter(std::string const &param, std::string_view const value) const
Definition: ConfigTree.cpp:151
Global vector based on Eigen vector.
Definition: EigenVector.h:25
std::size_t getID() const
Definition: Point3dWithID.h:62
A subset of nodes on a single mesh.
Definition: MeshSubset.h:26
bool isAxiallySymmetric() const
Definition: Mesh.h:133
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
Definition: Mesh.h:102
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition: Mesh.h:105
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
Definition: Mesh.h:84
std::size_t getID() const
Get id of the mesh.
Definition: Mesh.h:117
std::size_t getNumberOfNodes() const
Get the number of nodes.
Definition: Mesh.h:96
std::vector< Element const * > const & getElementsConnectedToNode(std::size_t node_id) const
Definition: Mesh.cpp:246
std::size_t getNumberOfElements() const
Get the number of elements.
Definition: Mesh.h:93
int getGlobalComponent(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 constexpr 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
std::unique_ptr< NumLib::LocalToGlobalIndexMap > _dof_table_boundary
Local dof table for the boundary mesh.
double interpolateToHigherOrderNode(GlobalVector const &x, int const var, int const comp, MeshLib::Node const &boundary_node) const
PythonBcData _bc_data
Auxiliary data used by the local assemblers.
std::vector< std::unique_ptr< PythonBoundaryConditionLocalAssemblerInterface > > _local_assemblers
Local assemblers for all elements of 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.
GlobalIndexType getDofIdx(std::size_t const boundary_node_id) const
void collectPrimaryVariables(std::vector< double > &primary_variables, MeshLib::Node const &boundary_node, GlobalVector const &x) const
PythonBoundaryCondition(PythonBcData &&bc_data, unsigned const integration_order, bool const flush_stdout, unsigned const bulk_mesh_dimension, NumLib::LocalToGlobalIndexMap const &dof_table_bulk)
void applyNaturalBC(const double t, std::vector< GlobalVector * > const &x, int const process_id, GlobalMatrix &K, GlobalVector &b, GlobalMatrix *Jac) override
void createLocalAssemblersPython(const unsigned dimension, std::vector< MeshLib::Element * > const &mesh_elements, NumLib::LocalToGlobalIndexMap const &dof_table, std::vector< std::unique_ptr< LocalAssemblerInterface > > &local_assemblers, NumLib::IntegrationOrder const integration_order, ExtraCtorArgs &&... extra_ctor_args)
std::unique_ptr< PythonBoundaryCondition > createPythonBoundaryCondition(BaseLib::ConfigTree const &config, MeshLib::Mesh const &boundary_mesh, NumLib::LocalToGlobalIndexMap const &dof_table_bulk, MeshLib::Mesh const &bulk_mesh, int const variable_id, int const component_id, unsigned const integration_order, unsigned const shapefunction_order, std::vector< std::reference_wrapper< ProcessVariable > > const &all_process_variables_for_this_process)
Creates a new PythonBoundaryCondition object.
void checkConsistency(NumLib::LocalToGlobalIndexMap const &dof_table, std::vector< std::reference_wrapper< ProcessLib::ProcessVariable > > const &pvs)
void initBCValues(NumLib::IndexValueVector< GlobalIndexType > &bc_values, std::size_t const nnodes)
std::vector< IndexType > ids
std::vector< double > values
static void executeMemberOnDereferenced(Method method, Container const &container, Args &&... args)
MeshLib::Mesh const & bc_or_st_mesh
The domain, where this BC or ST will be applied.
Definition: BcOrStData.h:46
std::vector< std::reference_wrapper< ProcessVariable > > const & all_process_variables_for_this_process
Definition: BcOrStData.h:51
BcOrStPythonSideInterface const *const bc_or_st_object
Python object computing BC or ST values.
Definition: BcOrStData.h:39