OGS
ConstraintDirichletBoundaryCondition.cpp
Go to the documentation of this file.
1
12
13#include <algorithm>
14#include <vector>
15
16#include "BaseLib/Logging.h"
17#include "MeshLib/MeshSearch/NodeSearch.h" // for getUniqueNodes
18#include "MeshLib/Node.h"
19#include "ParameterLib/Utils.h"
21
22namespace ProcessLib
23{
25 ParameterLib::Parameter<double> const& parameter,
26 NumLib::LocalToGlobalIndexMap const& dof_table_bulk, int const variable_id,
27 int const component_id, MeshLib::Mesh const& bc_mesh,
28 unsigned const integration_order, MeshLib::Mesh const& bulk_mesh,
29 double const constraint_threshold, bool const lower,
30 std::function<Eigen::Vector3d(std::size_t const, MathLib::Point3d const&,
31 double const,
32 std::vector<GlobalVector*> const&)>
33 getFlux)
34 : _parameter(parameter),
35 _variable_id(variable_id),
36 _component_id(component_id),
37 _bc_mesh(bc_mesh),
38 _integration_order(integration_order),
39 _constraint_threshold(constraint_threshold),
40 _lower(lower),
41 _bulk_mesh(bulk_mesh),
42 _getFlux(getFlux)
43{
44 if (variable_id >=
45 static_cast<int>(dof_table_bulk.getNumberOfVariables()) ||
46 component_id >=
47 dof_table_bulk.getNumberOfVariableComponents(variable_id))
48 {
50 "Variable id or component id too high. Actual values: ({:d}, "
51 "{:d}), maximum values: ({:d}, {:d}).",
52 variable_id, component_id, dof_table_bulk.getNumberOfVariables(),
53 dof_table_bulk.getNumberOfVariableComponents(variable_id));
54 }
55
56 std::vector<MeshLib::Node*> const& bc_nodes = _bc_mesh.getNodes();
57 DBUG(
58 "Found {:d} nodes for constraint Dirichlet BCs for the variable {:d} "
59 "and component {:d}",
60 bc_nodes.size(), variable_id, component_id);
61
62 MeshLib::MeshSubset bc_mesh_subset{_bc_mesh, bc_nodes};
63
64 // Create local DOF table from intersected mesh subsets for the given
65 // variable and component ids.
67 variable_id, {component_id}, std::move(bc_mesh_subset));
68
69 auto const& bc_elements(_bc_mesh.getElements());
70 _local_assemblers.resize(bc_elements.size());
71 _flux_values.resize(bc_elements.size());
72 // create _bulk_ids vector
73 auto const* bulk_element_ids = MeshLib::bulkElementIDs(_bc_mesh);
74 auto const* bulk_node_ids = MeshLib::bulkNodeIDs(_bc_mesh);
75 auto const& bulk_nodes = bulk_mesh.getNodes();
76
77 auto get_bulk_element_face_id =
78 [&](auto const bulk_element_id, MeshLib::Element const* bc_elem)
79 {
80 auto const* bulk_elem = _bulk_mesh.getElement(bulk_element_id);
81 std::array<MeshLib::Node const*, 3> nodes{
82 {bulk_nodes[(*bulk_node_ids)[bc_elem->getNode(0)->getID()]],
83 bulk_nodes[(*bulk_node_ids)[bc_elem->getNode(1)->getID()]],
84 bulk_nodes[(*bulk_node_ids)[bc_elem->getNode(2)->getID()]]}};
85 return bulk_elem->identifyFace(nodes.data());
86 };
87
88 _bulk_ids.reserve(bc_elements.size());
89 std::transform(
90 begin(bc_elements), end(bc_elements), std::back_inserter(_bulk_ids),
91 [&](auto const* bc_element)
92 {
93 auto const bulk_element_id =
94 (*bulk_element_ids)[bc_element->getID()];
95 return std::make_pair(
96 bulk_element_id,
97 get_bulk_element_face_id(bulk_element_id, bc_element));
98 });
99
100 const int shape_function_order = 1;
101
105 shape_function_order, _local_assemblers,
108}
109
111 double const t, std::vector<GlobalVector*> const& x,
112 int const /*process_id*/)
113{
114 DBUG(
115 "ConstraintDirichletBoundaryCondition::preTimestep: computing flux "
116 "constraints");
117 for (auto const id : _bc_mesh.getElements() | MeshLib::views::ids)
118 {
119 _flux_values[id] = _local_assemblers[id]->integrate(
120 x, t,
121 [this](std::size_t const element_id, MathLib::Point3d const& pnt,
122 double const t, std::vector<GlobalVector*> const& x)
123 { return _getFlux(element_id, pnt, t, x); });
124 }
125}
126
128 const double t, const GlobalVector& /*x*/,
130{
132
133 bc_values.ids.clear();
134 bc_values.values.clear();
135
136 std::vector<std::pair<GlobalIndexType, double>> tmp_bc_values;
137
138 auto isFlux = [&](const std::size_t element_id)
139 {
140 return _lower ? _flux_values[element_id] < _constraint_threshold
141 : _flux_values[element_id] > _constraint_threshold;
142 };
143
144 for (auto const* boundary_element : _bc_mesh.getElements())
145 {
146 // check if the boundary element is active
147 if (isFlux(boundary_element->getID()))
148 {
149 continue;
150 }
151
152 // loop over the boundary element nodes and set the dirichlet value
153 unsigned const number_nodes = boundary_element->getNumberOfNodes();
154 for (unsigned i = 0; i < number_nodes; ++i)
155 {
156 auto const id = boundary_element->getNode(i)->getID();
157 pos.setAll(id, boundary_element->getID(), {}, {});
158
160 id);
161 // TODO: that might be slow, but only done once
162 const auto g_idx = _dof_table_boundary->getGlobalIndex(
165 {
166 continue;
167 }
168 // For the DDC approach (e.g. with PETSc option), the negative
169 // index of g_idx means that the entry by that index is a ghost one,
170 // which should be dropped. Especially for PETSc routines
171 // MatZeroRows and MatZeroRowsColumns, which are called to apply the
172 // Dirichlet BC, the negative index is not accepted like other
173 // matrix or vector PETSc routines. Therefore, the following
174 // if-condition is applied.
175 if (g_idx >= 0)
176 {
177 tmp_bc_values.emplace_back(g_idx, _parameter(t, pos).front());
178 }
179 }
180 }
181
182 if (tmp_bc_values.empty())
183 {
184 DBUG("The domain on which the boundary is defined is empty");
185 return;
186 }
187
188 // store unique pairs of node id and value with respect to the node id in
189 // the bc_values vector. The values related to the same node id are
190 // averaged.
191 // first: sort the (node id, value) pairs according to the node id
192 std::sort(tmp_bc_values.begin(), tmp_bc_values.end(),
193 [](std::pair<GlobalIndexType, double> const& a,
194 std::pair<GlobalIndexType, double> const& b)
195 { return a.first < b.first; });
196 // second: average the values over equal node id ranges
197 unsigned cnt = 1;
198 GlobalIndexType current_id = tmp_bc_values.begin()->first;
199 double sum = tmp_bc_values.begin()->second;
200 for (auto const& tmp_bc_value : tmp_bc_values)
201 {
202 if (tmp_bc_value.first == current_id)
203 {
204 cnt++;
205 sum += tmp_bc_value.second;
206 }
207 else
208 {
209 bc_values.ids.emplace_back(current_id);
210 bc_values.values.emplace_back(sum / cnt);
211 cnt = 1;
212 sum = tmp_bc_value.second;
213 current_id = tmp_bc_value.first;
214 }
215 }
216 bc_values.ids.emplace_back(current_id);
217 bc_values.values.emplace_back(sum / cnt);
218
219 DBUG("Found {:d} constraint dirichlet boundary condition values.",
220 bc_values.ids.size());
221}
222
223std::unique_ptr<ConstraintDirichletBoundaryCondition>
225 BaseLib::ConfigTree const& config, MeshLib::Mesh const& bc_mesh,
226 NumLib::LocalToGlobalIndexMap const& dof_table_bulk, int const variable_id,
227 unsigned const integration_order, int const component_id,
228 std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
229 Process const& constraining_process)
230{
231 DBUG("Constructing ConstraintDirichletBoundaryCondition from config.");
233 config.checkConfigParameter("type", "ConstraintDirichlet");
234
235 auto const constraint_type =
237 config.getConfigParameter<std::string>("constraint_type");
238 if (constraint_type != "Flux")
239 {
240 OGS_FATAL("The constraint type is '{:s}', but has to be 'Flux'.",
241 constraint_type);
242 }
243
244 // Todo (TF) Open question: How to specify which getFlux function should be
245 // used for the constraint calculation?
246 auto const constraining_process_variable =
248 config.getConfigParameter<std::string>("constraining_process_variable");
249
250 if (!constraining_process.isMonolithicSchemeUsed())
251 {
252 OGS_FATAL(
253 "The constraint dirichlet boundary condition is implemented only "
254 "for monolithic implemented processes.");
255 }
256 const int process_id = 0;
257 auto process_variables =
258 constraining_process.getProcessVariables(process_id);
259 auto constraining_pv =
260 std::find_if(process_variables.cbegin(), process_variables.cend(),
261 [&constraining_process_variable](ProcessVariable const& pv)
262 { return pv.getName() == constraining_process_variable; });
263 if (constraining_pv == std::end(process_variables))
264 {
265 auto const& constraining_process_variable_name =
266 process_variables[variable_id].get().getName();
267 OGS_FATAL(
268 "<constraining_process_variable> in process variable name '{:s}' "
269 "at geometry 'TODO' : The constraining process variable is set as "
270 "'{:s}', but this is not specified in the project file.",
271 constraining_process_variable_name,
272 constraining_process_variable);
273 }
274
275 auto const constraint_threshold =
277 config.getConfigParameter<double>("constraint_threshold");
278
279 auto const constraint_direction_string =
281 config.getConfigParameter<std::string>("constraint_direction");
282 if (constraint_direction_string != "greater" &&
283 constraint_direction_string != "lower")
284 {
285 OGS_FATAL(
286 "The constraint direction is '{:s}', but has to be either "
287 "'greater' or 'lower'.",
288 constraint_direction_string);
289 }
290 bool const lower = constraint_direction_string == "lower";
291
293 auto const param_name = config.getConfigParameter<std::string>("parameter");
294 DBUG("Using parameter {:s}", param_name);
295
296 auto& param = ParameterLib::findParameter<double>(param_name, parameters, 1,
297 &bc_mesh);
298
299 // In case of partitioned mesh the boundary could be empty, i.e. there is no
300 // boundary condition.
301#ifdef USE_PETSC
302 // This can be extracted to createBoundaryCondition() but then the config
303 // parameters are not read and will cause an error.
304 // TODO (naumov): Add a function to ConfigTree for skipping the tags of the
305 // subtree and move the code up in createBoundaryCondition().
306 if (bc_mesh.getDimension() == 0 && bc_mesh.getNumberOfNodes() == 0 &&
307 bc_mesh.getNumberOfElements() == 0)
308 {
309 return nullptr;
310 }
311#endif // USE_PETSC
312
313 return std::make_unique<ConstraintDirichletBoundaryCondition>(
314 param, dof_table_bulk, variable_id, component_id, bc_mesh,
315 integration_order, constraining_process.getMesh(), constraint_threshold,
316 lower,
317 [&constraining_process](std::size_t const element_id,
318 MathLib::Point3d const& pnt, double const t,
319 std::vector<GlobalVector*> const& x)
320 { return constraining_process.getFlux(element_id, pnt, t, x); });
321}
322
323} // 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
Definition of the Node class.
T getConfigParameter(std::string const &param) const
void checkConfigParameter(std::string const &param, std::string_view const value) const
Global vector based on Eigen vector.
Definition EigenVector.h:25
virtual unsigned identifyFace(Node const *nodes[3]) const =0
Returns the ID of a face given an array of nodes.
A subset of nodes on a single mesh.
Definition MeshSubset.h:26
bool isAxiallySymmetric() const
Definition Mesh.h:137
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
Definition Mesh.h:106
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition Mesh.h:109
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
Definition Mesh.h:88
std::size_t getID() const
Get id of the mesh.
Definition Mesh.h:121
const Element * getElement(std::size_t idx) const
Get the element with the given index.
Definition Mesh.h:94
std::size_t getNumberOfNodes() const
Get the number of nodes.
Definition Mesh.h:100
std::size_t getNumberOfElements() const
Get the number of elements.
Definition Mesh.h:97
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
void setAll(std::optional< std::size_t > const &node_id, std::optional< std::size_t > const &element_id, std::optional< unsigned > const &integration_point, std::optional< MathLib::Point3d > const &coordinates)
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::function< Eigen::Vector3d(std::size_t const, MathLib::Point3d const &, double const, std::vector< GlobalVector * > const &)> _getFlux
The function _getFlux calculates the flux through the boundary element.
std::vector< std::unique_ptr< ConstraintDirichletBoundaryConditionLocalAssemblerInterface > > _local_assemblers
Local assemblers for each boundary element.
std::unique_ptr< NumLib::LocalToGlobalIndexMap > _dof_table_boundary
void preTimestep(double const t, std::vector< GlobalVector * > const &x, int const process_id) override
std::vector< std::pair< std::size_t, unsigned > > _bulk_ids
std::vector< double > _flux_values
Stores the results of the flux computations per boundary element.
ConstraintDirichletBoundaryCondition(ParameterLib::Parameter< double > const &parameter, NumLib::LocalToGlobalIndexMap const &dof_table_bulk, int const variable_id, int const component_id, MeshLib::Mesh const &bc_mesh, unsigned const integration_order, MeshLib::Mesh const &bulk_mesh, double const constraint_threshold, bool const lower, std::function< Eigen::Vector3d(std::size_t const, MathLib::Point3d const &, double const, std::vector< GlobalVector * > const &)> getFlux)
unsigned const _integration_order
Integration order for integration over the lower-dimensional elements.
MeshLib::Mesh & getMesh() const
Definition Process.h:153
virtual bool isMonolithicSchemeUsed() const
Definition Process.h:102
std::vector< std::vector< std::reference_wrapper< ProcessVariable > > > const & getProcessVariables() const
Definition Process.h:156
constexpr ranges::views::view_closure ids
For an element of a range view return its id.
Definition Mesh.h:225
PropertyVector< std::size_t > const * bulkElementIDs(Mesh const &mesh)
Definition Mesh.cpp:300
PropertyVector< std::size_t > const * bulkNodeIDs(Mesh const &mesh)
Definition Mesh.cpp:292
OGS_NO_DANGLING Parameter< ParameterDataType > & findParameter(std::string const &parameter_name, std::vector< std::unique_ptr< ParameterBase > > const &parameters, int const num_components, MeshLib::Mesh const *const mesh=nullptr)
Definition Utils.h:102
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, NumLib::IntegrationOrder const integration_order, ExtraCtorArgs &&... extra_ctor_args)
std::unique_ptr< ConstraintDirichletBoundaryCondition > createConstraintDirichletBoundaryCondition(BaseLib::ConfigTree const &config, MeshLib::Mesh const &bc_mesh, NumLib::LocalToGlobalIndexMap const &dof_table_bulk, int const variable_id, unsigned const integration_order, int const component_id, std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const &parameters, Process const &constraining_process)
std::vector< IndexType > ids
std::vector< double > values