28 unsigned const integration_order,
MeshLib::Mesh const& bulk_mesh,
29 double const constraint_threshold,
bool const lower,
32 std::vector<GlobalVector*>
const&)>
34 : _parameter(parameter),
35 _variable_id(variable_id),
36 _component_id(component_id),
38 _integration_order(integration_order),
39 _constraint_threshold(constraint_threshold),
41 _bulk_mesh(bulk_mesh),
50 "Variable id or component id too high. Actual values: ({:d}, "
51 "{:d}), maximum values: ({:d}, {:d}).",
58 "Found {:d} nodes for constraint Dirichlet BCs for the variable {:d} "
60 bc_nodes.size(), variable_id, component_id);
67 variable_id, {component_id}, std::move(bc_mesh_subset)));
73 auto const* bulk_element_ids =
76 auto const* bulk_node_ids =
79 auto const& bulk_nodes = bulk_mesh.
getNodes();
81 auto get_bulk_element_face_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()]]}};
94 begin(bc_elements), end(bc_elements), std::back_inserter(
_bulk_ids),
95 [&](
auto const* bc_element)
97 auto const bulk_element_id =
98 (*bulk_element_ids)[bc_element->getID()];
99 return std::make_pair(
101 get_bulk_element_face_id(bulk_element_id, bc_element));
104 const int shape_function_order = 1;
114 double const t, std::vector<GlobalVector*>
const& x,
118 "ConstraintDirichletBoundaryCondition::preTimestep: computing flux "
125 [
this](std::size_t
const element_id,
127 std::vector<GlobalVector*>
const& x)
128 {
return _getFlux(element_id, pnt, t, x); });
138 bc_values.
ids.clear();
141 std::vector<std::pair<GlobalIndexType, double>> tmp_bc_values;
143 auto isFlux = [&](
const std::size_t element_id)
152 if (isFlux(boundary_element->getID()))
158 unsigned const number_nodes = boundary_element->getNumberOfNodes();
159 for (
unsigned i = 0; i < number_nodes; ++i)
161 auto const id = boundary_element->getNode(i)->getID();
162 pos.
setAll(
id, boundary_element->getID(), {}, {});
182 tmp_bc_values.emplace_back(g_idx,
_parameter(t, pos).front());
187 if (tmp_bc_values.empty())
189 DBUG(
"The domain on which the boundary is defined is empty");
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; });
204 double sum = tmp_bc_values.begin()->second;
205 for (
auto const& tmp_bc_value : tmp_bc_values)
207 if (tmp_bc_value.first == current_id)
210 sum += tmp_bc_value.second;
214 bc_values.
ids.emplace_back(current_id);
215 bc_values.
values.emplace_back(sum / cnt);
217 sum = tmp_bc_value.second;
218 current_id = tmp_bc_value.first;
221 bc_values.
ids.emplace_back(current_id);
222 bc_values.
values.emplace_back(sum / cnt);
224 DBUG(
"Found {:d} constraint dirichlet boundary condition values.",
225 bc_values.
ids.size());
228 std::unique_ptr<ConstraintDirichletBoundaryCondition>
232 unsigned const integration_order,
int const component_id,
233 std::vector<std::unique_ptr<ParameterLib::ParameterBase>>
const& parameters,
234 Process const& constraining_process)
236 DBUG(
"Constructing ConstraintDirichletBoundaryCondition from config.");
240 auto const constraint_type =
243 if (constraint_type !=
"Flux")
245 OGS_FATAL(
"The constraint type is '{:s}', but has to be 'Flux'.",
251 auto const constraining_process_variable =
258 "The constraint dirichlet boundary condition is implemented only "
259 "for monolithic implemented processes.");
261 const int process_id = 0;
262 auto process_variables =
264 auto constraining_pv =
265 std::find_if(process_variables.cbegin(), process_variables.cend(),
267 { return pv.getName() == constraining_process_variable; });
268 if (constraining_pv == std::end(process_variables))
270 auto const& constraining_process_variable_name =
271 process_variables[variable_id].get().getName();
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);
280 auto const constraint_threshold =
284 auto const constraint_direction_string =
287 if (constraint_direction_string !=
"greater" &&
288 constraint_direction_string !=
"lower")
291 "The constraint direction is '{:s}', but has to be either "
292 "'greater' or 'lower'.",
293 constraint_direction_string);
295 bool const lower = constraint_direction_string ==
"lower";
299 DBUG(
"Using parameter {:s}", param_name);
301 auto& param = ParameterLib::findParameter<double>(param_name, parameters, 1,
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,
322 [&constraining_process](std::size_t
const element_id,
324 std::vector<GlobalVector*>
const& x)
325 { return constraining_process.getFlux(element_id, pnt, t, x); });
GlobalMatrix::IndexType GlobalIndexType
void DBUG(char const *fmt, Args const &... args)
Definition of the Node class.
void checkConfigParameter(std::string const ¶m, T const &value) const
T getConfigParameter(std::string const ¶m) const
Global vector based on Eigen vector.
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.
bool isAxiallySymmetric() const
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
std::size_t getID() const
Get id of the mesh.
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Properties & getProperties()
std::size_t getNumberOfNodes() const
Get the number of nodes.
const Element * getElement(std::size_t idx) const
Get the element with the given index.
std::size_t getNumberOfElements() const
Get the number of elements.
PropertyVector< T > const * getPropertyVector(std::string const &name) const
int getNumberOfVariables() const
int getNumberOfVariableComponents(int variable_id) const
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)
ParameterLib::Parameter< double > const & _parameter
void getEssentialBCValues(const double t, const GlobalVector &x, NumLib::IndexValueVector< GlobalIndexType > &bc_values) const override
Writes the values of essential BCs to bc_values.
double const _constraint_threshold
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
MeshLib::Mesh const & _bc_mesh
MeshLib::Mesh const & _bulk_mesh
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 ¶meter, 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
MeshLib::Mesh & getMesh() const
bool isMonolithicSchemeUsed() const
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 ¶meters, Process const &constraining_process)
std::vector< IndexType > ids
std::vector< double > values