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 
22 namespace 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  {
49  OGS_FATAL(
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 =
75  "bulk_element_ids", MeshLib::MeshItemType::Cell, 1);
76  auto const* bulk_node_ids =
78  "bulk_node_ids", MeshLib::MeshItemType::Node, 1);
79  auto const& bulk_nodes = bulk_mesh.getNodes();
80 
81  auto get_bulk_element_face_id =
82  [&](auto const bulk_element_id, MeshLib::Element const* bc_elem)
83  {
84  auto const* bulk_elem = _bulk_mesh.getElement(bulk_element_id);
85  std::array<MeshLib::Node const*, 3> nodes{
86  {bulk_nodes[(*bulk_node_ids)[bc_elem->getNode(0)->getID()]],
87  bulk_nodes[(*bulk_node_ids)[bc_elem->getNode(1)->getID()]],
88  bulk_nodes[(*bulk_node_ids)[bc_elem->getNode(2)->getID()]]}};
89  return bulk_elem->identifyFace(nodes.data());
90  };
91 
92  _bulk_ids.reserve(bc_elements.size());
93  std::transform(
94  begin(bc_elements), end(bc_elements), std::back_inserter(_bulk_ids),
95  [&](auto const* bc_element)
96  {
97  auto const bulk_element_id =
98  (*bulk_element_ids)[bc_element->getID()];
99  return std::make_pair(
100  bulk_element_id,
101  get_bulk_element_face_id(bulk_element_id, bc_element));
102  });
103 
104  const int shape_function_order = 1;
105 
109  shape_function_order, _local_assemblers, _bc_mesh.isAxiallySymmetric(),
110  _integration_order, bulk_mesh, _bulk_ids);
111 }
112 
114  double const t, std::vector<GlobalVector*> const& x,
115  int const /*process_id*/)
116 {
117  DBUG(
118  "ConstraintDirichletBoundaryCondition::preTimestep: computing flux "
119  "constraints");
120  for (auto const* boundary_element : _bc_mesh.getElements())
121  {
122  _flux_values[boundary_element->getID()] =
123  _local_assemblers[boundary_element->getID()]->integrate(
124  x, t,
125  [this](std::size_t const element_id,
126  MathLib::Point3d const& pnt, double const t,
127  std::vector<GlobalVector*> const& x)
128  { return _getFlux(element_id, pnt, t, x); });
129  }
130 }
131 
133  const double t, const GlobalVector& /*x*/,
135 {
137 
138  bc_values.ids.clear();
139  bc_values.values.clear();
140 
141  std::vector<std::pair<GlobalIndexType, double>> tmp_bc_values;
142 
143  auto isFlux = [&](const std::size_t element_id)
144  {
145  return _lower ? _flux_values[element_id] < _constraint_threshold
146  : _flux_values[element_id] > _constraint_threshold;
147  };
148 
149  for (auto const* boundary_element : _bc_mesh.getElements())
150  {
151  // check if the boundary element is active
152  if (isFlux(boundary_element->getID()))
153  {
154  continue;
155  }
156 
157  // loop over the boundary element nodes and set the dirichlet value
158  unsigned const number_nodes = boundary_element->getNumberOfNodes();
159  for (unsigned i = 0; i < number_nodes; ++i)
160  {
161  auto const id = boundary_element->getNode(i)->getID();
162  pos.setAll(id, boundary_element->getID(), {}, {});
163 
165  id);
166  // TODO: that might be slow, but only done once
167  const auto g_idx = _dof_table_boundary->getGlobalIndex(
169  if (g_idx == NumLib::MeshComponentMap::nop)
170  {
171  continue;
172  }
173  // For the DDC approach (e.g. with PETSc option), the negative
174  // index of g_idx means that the entry by that index is a ghost one,
175  // which should be dropped. Especially for PETSc routines
176  // MatZeroRows and MatZeroRowsColumns, which are called to apply the
177  // Dirichlet BC, the negative index is not accepted like other
178  // matrix or vector PETSc routines. Therefore, the following
179  // if-condition is applied.
180  if (g_idx >= 0)
181  {
182  tmp_bc_values.emplace_back(g_idx, _parameter(t, pos).front());
183  }
184  }
185  }
186 
187  if (tmp_bc_values.empty())
188  {
189  DBUG("The domain on which the boundary is defined is empty");
190  return;
191  }
192 
193  // store unique pairs of node id and value with respect to the node id in
194  // the bc_values vector. The values related to the same node id are
195  // averaged.
196  // first: sort the (node id, value) pairs according to the node id
197  std::sort(tmp_bc_values.begin(), tmp_bc_values.end(),
198  [](std::pair<GlobalIndexType, double> const& a,
199  std::pair<GlobalIndexType, double> const& b)
200  { return a.first < b.first; });
201  // second: average the values over equal node id ranges
202  unsigned cnt = 1;
203  GlobalIndexType current_id = tmp_bc_values.begin()->first;
204  double sum = tmp_bc_values.begin()->second;
205  for (auto const& tmp_bc_value : tmp_bc_values)
206  {
207  if (tmp_bc_value.first == current_id)
208  {
209  cnt++;
210  sum += tmp_bc_value.second;
211  }
212  else
213  {
214  bc_values.ids.emplace_back(current_id);
215  bc_values.values.emplace_back(sum / cnt);
216  cnt = 1;
217  sum = tmp_bc_value.second;
218  current_id = tmp_bc_value.first;
219  }
220  }
221  bc_values.ids.emplace_back(current_id);
222  bc_values.values.emplace_back(sum / cnt);
223 
224  DBUG("Found {:d} constraint dirichlet boundary condition values.",
225  bc_values.ids.size());
226 }
227 
228 std::unique_ptr<ConstraintDirichletBoundaryCondition>
230  BaseLib::ConfigTree const& config, MeshLib::Mesh const& bc_mesh,
231  NumLib::LocalToGlobalIndexMap const& dof_table_bulk, int const variable_id,
232  unsigned const integration_order, int const component_id,
233  std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
234  Process const& constraining_process)
235 {
236  DBUG("Constructing ConstraintDirichletBoundaryCondition from config.");
238  config.checkConfigParameter("type", "ConstraintDirichlet");
239 
240  auto const constraint_type =
242  config.getConfigParameter<std::string>("constraint_type");
243  if (constraint_type != "Flux")
244  {
245  OGS_FATAL("The constraint type is '{:s}', but has to be 'Flux'.",
246  constraint_type);
247  }
248 
249  // Todo (TF) Open question: How to specify which getFlux function should be
250  // used for the constraint calculation?
251  auto const constraining_process_variable =
253  config.getConfigParameter<std::string>("constraining_process_variable");
254 
255  if (!constraining_process.isMonolithicSchemeUsed())
256  {
257  OGS_FATAL(
258  "The constraint dirichlet boundary condition is implemented only "
259  "for monolithic implemented processes.");
260  }
261  const int process_id = 0;
262  auto process_variables =
263  constraining_process.getProcessVariables(process_id);
264  auto constraining_pv =
265  std::find_if(process_variables.cbegin(), process_variables.cend(),
266  [&constraining_process_variable](ProcessVariable const& pv)
267  { return pv.getName() == constraining_process_variable; });
268  if (constraining_pv == std::end(process_variables))
269  {
270  auto const& constraining_process_variable_name =
271  process_variables[variable_id].get().getName();
272  OGS_FATAL(
273  "<constraining_process_variable> in process variable name '{:s}' "
274  "at geometry 'TODO' : The constraining process variable is set as "
275  "'{:s}', but this is not specified in the project file.",
276  constraining_process_variable_name,
277  constraining_process_variable);
278  }
279 
280  auto const constraint_threshold =
282  config.getConfigParameter<double>("constraint_threshold");
283 
284  auto const constraint_direction_string =
286  config.getConfigParameter<std::string>("constraint_direction");
287  if (constraint_direction_string != "greater" &&
288  constraint_direction_string != "lower")
289  {
290  OGS_FATAL(
291  "The constraint direction is '{:s}', but has to be either "
292  "'greater' or 'lower'.",
293  constraint_direction_string);
294  }
295  bool const lower = constraint_direction_string == "lower";
296 
298  auto const param_name = config.getConfigParameter<std::string>("parameter");
299  DBUG("Using parameter {:s}", param_name);
300 
301  auto& param = ParameterLib::findParameter<double>(param_name, parameters, 1,
302  &bc_mesh);
303 
304  // In case of partitioned mesh the boundary could be empty, i.e. there is no
305  // boundary condition.
306 #ifdef USE_PETSC
307  // This can be extracted to createBoundaryCondition() but then the config
308  // parameters are not read and will cause an error.
309  // TODO (naumov): Add a function to ConfigTree for skipping the tags of the
310  // subtree and move the code up in createBoundaryCondition().
311  if (bc_mesh.getDimension() == 0 && bc_mesh.getNumberOfNodes() == 0 &&
312  bc_mesh.getNumberOfElements() == 0)
313  {
314  return nullptr;
315  }
316 #endif // USE_PETSC
317 
318  return std::make_unique<ConstraintDirichletBoundaryCondition>(
319  param, dof_table_bulk, variable_id, component_id, bc_mesh,
320  integration_order, constraining_process.getMesh(), constraint_threshold,
321  lower,
322  [&constraining_process](std::size_t const element_id,
323  MathLib::Point3d const& pnt, double const t,
324  std::vector<GlobalVector*> const& x)
325  { return constraining_process.getFlux(element_id, pnt, t, x); });
326 }
327 
328 } // namespace ProcessLib
#define OGS_FATAL(...)
Definition: Error.h:26
GlobalMatrix::IndexType GlobalIndexType
void DBUG(char const *fmt, Args const &... args)
Definition: Logging.h:27
Definition of the Node class.
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
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: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::size_t getID() const
Get id of the mesh.
Definition: Mesh.h:110
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
const Element * getElement(std::size_t idx) const
Get the element with the given index.
Definition: Mesh.h:77
std::size_t getNumberOfElements() const
Get the number of elements.
Definition: Mesh.h:86
PropertyVector< T > const * getPropertyVector(std::string const &name) 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
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::TemplatePoint< double, 3 >> 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.
std::vector< std::reference_wrapper< ProcessVariable > > const & getProcessVariables(const int process_id) const
Definition: Process.h:144
MeshLib::Mesh & getMesh() const
Definition: Process.h:142
bool isMonolithicSchemeUsed() const
Definition: Process.h:102
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< 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