OGS
ProcessLib::ConstraintDirichletBoundaryCondition Class Referencefinal

Detailed Description

The ConstraintDirichletBoundaryCondition class describes a Dirichlet-type boundary condition that is constant in space and time where the domain can shrink and grow within the simulation. The expected parameter in the passed configuration is "value" which, when not present defaults to zero.

Definition at line 25 of file ConstraintDirichletBoundaryCondition.h.

#include <ConstraintDirichletBoundaryCondition.h>

Inheritance diagram for ProcessLib::ConstraintDirichletBoundaryCondition:
[legend]
Collaboration diagram for ProcessLib::ConstraintDirichletBoundaryCondition:
[legend]

Public Member Functions

 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)
 
void preTimestep (double const t, std::vector< GlobalVector * > const &x, int const process_id) override
 
void getEssentialBCValues (const double t, const GlobalVector &x, NumLib::IndexValueVector< GlobalIndexType > &bc_values) const override
 Writes the values of essential BCs to bc_values. More...
 
- Public Member Functions inherited from ProcessLib::BoundaryCondition
virtual void applyNaturalBC (const double, std::vector< GlobalVector * > const &, int const, GlobalMatrix &, GlobalVector &, GlobalMatrix *)
 
virtual void postTimestep (const double, std::vector< GlobalVector * > const &, int const)
 
virtual ~BoundaryCondition ()=default
 

Private Attributes

ParameterLib::Parameter< double > const & _parameter
 
std::unique_ptr< NumLib::LocalToGlobalIndexMap_dof_table_boundary
 
int const _variable_id
 
int const _component_id
 
MeshLib::Mesh const & _bc_mesh
 
unsigned const _integration_order
 Integration order for integration over the lower-dimensional elements. More...
 
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. More...
 
std::vector< std::unique_ptr< ConstraintDirichletBoundaryConditionLocalAssemblerInterface > > _local_assemblers
 Local assemblers for each boundary element. More...
 
double const _constraint_threshold
 
bool const _lower
 
MeshLib::Mesh const & _bulk_mesh
 
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. More...
 

Constructor & Destructor Documentation

◆ ConstraintDirichletBoundaryCondition()

ProcessLib::ConstraintDirichletBoundaryCondition::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 
)
Parameters
parameterUsed for setting the values for the boundary condition.
dof_table_bulkThe bulk local to global index map is used to derive the local to global index map for the boundary.
variable_idThe variable id is needed to determine the global index.
component_idThe component id is needed to determine the global index.
bc_meshLower dimensional mesh the boundary condition is defined on. The bc_mesh must have the two PropertyVector objects 'bulk_element_ids' and 'bulk_node_ids' containing the corresponding information.
integration_orderOrder the order of integration used to compute the constraint.
bulk_meshThe FE mesh for the simulation.
constraint_thresholdThe threshold value used for the switch off/on decision.
lowerBoolean value used for the calculation of the constraint criterion, i.e., if lower is set to true the criterion 'calculated_value < constraint_threshold' is evaluated to switch on/off the boundary condition, else 'calculated_value > constraint_threshold' is evaluated.
getFluxThe function used for the flux calculation.
Note
The function has to be stored by value, else the process value is not captured properly.

Definition at line 24 of file ConstraintDirichletBoundaryCondition.cpp.

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.
66  _dof_table_boundary.reset(dof_table_bulk.deriveBoundaryConstrainedMap(
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 
107  ConstraintDirichletBoundaryConditionLocalAssembler>(
109  shape_function_order, _local_assemblers, _bc_mesh.isAxiallySymmetric(),
110  _integration_order, bulk_mesh, _bulk_ids);
111 }
#define OGS_FATAL(...)
Definition: Error.h:26
void DBUG(char const *fmt, Args const &... args)
Definition: Logging.h:27
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:27
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
const Element * getElement(std::size_t idx) const
Get the element with the given index.
Definition: Mesh.h:77
PropertyVector< T > const * getPropertyVector(std::string const &name) const
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
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.
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)

References _bc_mesh, _bulk_ids, _bulk_mesh, _dof_table_boundary, _flux_values, _integration_order, _local_assemblers, MeshLib::Cell, ProcessLib::BoundaryConditionAndSourceTerm::createLocalAssemblers(), DBUG(), NumLib::LocalToGlobalIndexMap::deriveBoundaryConstrainedMap(), MeshLib::Mesh::getDimension(), MeshLib::Mesh::getElement(), MeshLib::Mesh::getElements(), MeshLib::Mesh::getNodes(), NumLib::LocalToGlobalIndexMap::getNumberOfVariableComponents(), NumLib::LocalToGlobalIndexMap::getNumberOfVariables(), MeshLib::Mesh::getProperties(), MeshLib::Properties::getPropertyVector(), MeshLib::Element::identifyFace(), MeshLib::Mesh::isAxiallySymmetric(), MeshLib::Node, and OGS_FATAL.

Member Function Documentation

◆ getEssentialBCValues()

void ProcessLib::ConstraintDirichletBoundaryCondition::getEssentialBCValues ( const double  ,
const GlobalVector ,
NumLib::IndexValueVector< GlobalIndexType > &   
) const
overridevirtual

Writes the values of essential BCs to bc_values.

Reimplemented from ProcessLib::BoundaryCondition.

Definition at line 132 of file ConstraintDirichletBoundaryCondition.cpp.

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 }
GlobalMatrix::IndexType GlobalIndexType
std::size_t getID() const
Get id of the mesh.
Definition: Mesh.h:110
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)

References _bc_mesh, _component_id, _constraint_threshold, _dof_table_boundary, _flux_values, _lower, _parameter, _variable_id, DBUG(), MeshLib::Mesh::getElements(), MeshLib::Mesh::getID(), NumLib::IndexValueVector< IndexType >::ids, MeshLib::Node, NumLib::MeshComponentMap::nop, ParameterLib::SpatialPosition::setAll(), and NumLib::IndexValueVector< IndexType >::values.

◆ preTimestep()

void ProcessLib::ConstraintDirichletBoundaryCondition::preTimestep ( double const  t,
std::vector< GlobalVector * > const &  x,
int const  process_id 
)
overridevirtual

Reimplemented from ProcessLib::BoundaryCondition.

Definition at line 113 of file ConstraintDirichletBoundaryCondition.cpp.

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 }

References _bc_mesh, _flux_values, _getFlux, _local_assemblers, DBUG(), and MeshLib::Mesh::getElements().

Member Data Documentation

◆ _bc_mesh

MeshLib::Mesh const& ProcessLib::ConstraintDirichletBoundaryCondition::_bc_mesh
private

Vector of (lower-dimensional) boundary elements on which the boundary condition is defined.

Definition at line 82 of file ConstraintDirichletBoundaryCondition.h.

Referenced by ConstraintDirichletBoundaryCondition(), getEssentialBCValues(), and preTimestep().

◆ _bulk_ids

std::vector<std::pair<std::size_t, unsigned> > ProcessLib::ConstraintDirichletBoundaryCondition::_bulk_ids
private

The first item of the pair is the element id in the bulk mesh, the second item is the face id of the bulk element that is part of the boundary

Definition at line 90 of file ConstraintDirichletBoundaryCondition.h.

Referenced by ConstraintDirichletBoundaryCondition().

◆ _bulk_mesh

MeshLib::Mesh const& ProcessLib::ConstraintDirichletBoundaryCondition::_bulk_mesh
private

The mesh _bulk_mesh is the discretized domain the process(es) are defined on. It is needed to get values for the constraint calculation.

Definition at line 112 of file ConstraintDirichletBoundaryCondition.h.

Referenced by ConstraintDirichletBoundaryCondition().

◆ _component_id

int const ProcessLib::ConstraintDirichletBoundaryCondition::_component_id
private

Definition at line 78 of file ConstraintDirichletBoundaryCondition.h.

Referenced by getEssentialBCValues().

◆ _constraint_threshold

double const ProcessLib::ConstraintDirichletBoundaryCondition::_constraint_threshold
private

The threshold value used to the switch off/on the Dirichlet-type boundary condition.

Definition at line 102 of file ConstraintDirichletBoundaryCondition.h.

Referenced by getEssentialBCValues().

◆ _dof_table_boundary

std::unique_ptr<NumLib::LocalToGlobalIndexMap> ProcessLib::ConstraintDirichletBoundaryCondition::_dof_table_boundary
private

Local dof table, a subset of the global one restricted to the participating number of elements of the boundary condition.

Definition at line 75 of file ConstraintDirichletBoundaryCondition.h.

Referenced by ConstraintDirichletBoundaryCondition(), and getEssentialBCValues().

◆ _flux_values

std::vector<double> ProcessLib::ConstraintDirichletBoundaryCondition::_flux_values
private

Stores the results of the flux computations per boundary element.

Definition at line 93 of file ConstraintDirichletBoundaryCondition.h.

Referenced by ConstraintDirichletBoundaryCondition(), getEssentialBCValues(), and preTimestep().

◆ _getFlux

std::function<Eigen::Vector3d(std::size_t const, MathLib::Point3d const&, double const, std::vector<GlobalVector*> const&)> ProcessLib::ConstraintDirichletBoundaryCondition::_getFlux
private

The function _getFlux calculates the flux through the boundary element.

Definition at line 118 of file ConstraintDirichletBoundaryCondition.h.

Referenced by preTimestep().

◆ _integration_order

unsigned const ProcessLib::ConstraintDirichletBoundaryCondition::_integration_order
private

Integration order for integration over the lower-dimensional elements.

Definition at line 85 of file ConstraintDirichletBoundaryCondition.h.

Referenced by ConstraintDirichletBoundaryCondition().

◆ _local_assemblers

std::vector<std::unique_ptr< ConstraintDirichletBoundaryConditionLocalAssemblerInterface> > ProcessLib::ConstraintDirichletBoundaryCondition::_local_assemblers
private

Local assemblers for each boundary element.

Definition at line 98 of file ConstraintDirichletBoundaryCondition.h.

Referenced by ConstraintDirichletBoundaryCondition(), and preTimestep().

◆ _lower

bool const ProcessLib::ConstraintDirichletBoundaryCondition::_lower
private

The boolean value lower is used for the calculation of the constraint criterion, i.e., if lower is set to true the criterion 'calculated_value < constraint_threshold' is evaluated to switch on/off the boundary condition, else 'calculated_value > constraint_threshold' is evaluated.

Definition at line 108 of file ConstraintDirichletBoundaryCondition.h.

Referenced by getEssentialBCValues().

◆ _parameter

ParameterLib::Parameter<double> const& ProcessLib::ConstraintDirichletBoundaryCondition::_parameter
private

Definition at line 71 of file ConstraintDirichletBoundaryCondition.h.

Referenced by getEssentialBCValues().

◆ _variable_id

int const ProcessLib::ConstraintDirichletBoundaryCondition::_variable_id
private

Definition at line 77 of file ConstraintDirichletBoundaryCondition.h.

Referenced by getEssentialBCValues().


The documentation for this class was generated from the following files: