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