diff --git a/Applications/CLI/Tests.cmake b/Applications/CLI/Tests.cmake index 78616519c26f312b370a7ff7ec3d1d220d17a6be..8a912e19c02ac34f5bebf76469127914d07d02b3 100644 --- a/Applications/CLI/Tests.cmake +++ b/Applications/CLI/Tests.cmake @@ -134,6 +134,30 @@ if(NOT OGS_USE_MPI) DIFF_DATA line_1_line_${mesh_size}.vtu line_${mesh_size}_neumann_pcs_0_ts_1_t_1.000000.vtu D1_left_N1_right pressure ) + + AddTest( + NAME GroundWaterFlowProcess_line_1_Robin_Right_Picard_${mesh_size} + PATH Elliptic/line_1_GroundWaterFlow + EXECUTABLE ogs + EXECUTABLE_ARGS line_${mesh_size}_robin_right_picard.prj + WRAPPER time + TESTER vtkdiff + ABSTOL 1e-14 RELTOL 1e-14 + DIFF_DATA + line_1_line_${mesh_size}.vtu line_${mesh_size}_robin_right_picard_pcs_0_ts_1_t_1.000000.vtu D1_left_N1_right pressure + ) + + AddTest( + NAME GroundWaterFlowProcess_line_1_Robin_Left_Picard_${mesh_size} + PATH Elliptic/line_1_GroundWaterFlow + EXECUTABLE ogs + EXECUTABLE_ARGS line_${mesh_size}_robin_left_picard.prj + WRAPPER time + TESTER vtkdiff + ABSTOL 1e-14 RELTOL 1e-14 + DIFF_DATA + line_1_line_${mesh_size}.vtu line_${mesh_size}_robin_left_picard_pcs_0_ts_1_t_1.000000.vtu D1_left_N1_right pressure + ) endforeach() # Some Neumann BC tests diff --git a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp index 1678292f7de154937b7c04deb08eff5a7b1f2ce9..65ce716963ab88aa94e01528dec4966cacfc9777 100644 --- a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp +++ b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp @@ -13,6 +13,23 @@ #include "BoundaryConditionConfig.h" #include "UniformDirichletBoundaryCondition.h" #include "UniformNeumannBoundaryCondition.h" +#include "UniformRobinBoundaryCondition.h" + +static std::vector<MeshLib::Element*> getClonedElements( + MeshGeoToolsLib::BoundaryElementsSearcher& boundary_element_searcher, + GeoLib::GeoObject const& geometry) +{ + std::vector<MeshLib::Element*> elements = + boundary_element_searcher.getBoundaryElements(geometry); + + // Deep copy all the elements, because the searcher might destroy the + // originals. Store pointers to the copies in the elements vector (i.e., + // in-place modification). + for (auto& e : elements) + e = e->clone(); + + return elements; +} namespace ProcessLib { @@ -44,18 +61,18 @@ std::unique_ptr<BoundaryCondition> createBoundaryCondition( } else if (type == "UniformNeumann") { - std::vector<MeshLib::Element*> elements = - boundary_element_searcher.getBoundaryElements(config.geometry); - - // Deep copy all the elements, because the searcher might destroy the - // originals. Store pointers to the copies in the elements vector (i.e., - // in-place modification). - for (auto& e : elements) - e = e->clone(); - return createUniformNeumannBoundaryCondition( - config.config, std::move(elements), dof_table, variable_id, - config.component_id, integration_order, mesh.getDimension()); + config.config, + getClonedElements(boundary_element_searcher, config.geometry), + dof_table, variable_id, config.component_id, integration_order, + mesh.getDimension()); + } + else if (type == "UniformRobin") { + return createUniformRobinBoundaryCondition( + config.config, + getClonedElements(boundary_element_searcher, config.geometry), + dof_table, variable_id, config.component_id, integration_order, + mesh.getDimension()); } else { diff --git a/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.cpp new file mode 100644 index 0000000000000000000000000000000000000000..95bc0c90ca50bc7c9661a73fdb207e374efbe7e4 --- /dev/null +++ b/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.cpp @@ -0,0 +1,38 @@ +/** + * \copyright + * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#include "UniformRobinBoundaryCondition.h" + +namespace ProcessLib +{ +std::unique_ptr<UniformRobinBoundaryCondition> +createUniformRobinBoundaryCondition( + BaseLib::ConfigTree const& config, + std::vector<MeshLib::Element*>&& elements, + NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id, + int const component_id, unsigned const integration_order, + unsigned const global_dim) +{ + DBUG("Constructing RobinBcConfig from config."); + //! \ogs_file_param{boundary_condition__type} + config.checkConfigParameter("type", "UniformRobin"); + + //! \ogs_file_param{boundary_condition__UniformRobin__alpha} + double const alpha = config.getConfigParameter<double>("alpha"); + //! \ogs_file_param{boundary_condition__UniformRobin__u_0} + double const u_0 = config.getConfigParameter<double>("u_0"); + + return std::unique_ptr<UniformRobinBoundaryCondition>( + new UniformRobinBoundaryCondition( + integration_order, dof_table, variable_id, component_id, global_dim, + std::move(elements), + UniformRobinBoundaryConditionData{alpha, u_0})); +} + +} // ProcessLib diff --git a/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.h b/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.h new file mode 100644 index 0000000000000000000000000000000000000000..87ef77b8825768b1bfcc500520c209db921940c2 --- /dev/null +++ b/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.h @@ -0,0 +1,48 @@ +/** + * \copyright + * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#ifndef PROCESSLIB_UNIFORMROBINBOUNDARYCONDITION_H +#define PROCESSLIB_UNIFORMROBINBOUNDARYCONDITION_H + +#include "GenericNaturalBoundaryCondition.h" +#include "UniformRobinBoundaryConditionLocalAssembler.h" + +namespace MeshGeoToolsLib +{ +class BoundaryElementsSearcher; +} + +namespace ProcessLib +{ +using UniformRobinBoundaryCondition = GenericNaturalBoundaryCondition< + UniformRobinBoundaryConditionData, + UniformRobinBoundaryConditionLocalAssembler>; + +/*! Creates a new uniform Robin boundary condition from the given data. + * + * The Robin boundary condition is given in the form + * \f$ \alpha \cdot [ u_0 - u(x) ] \f$, + * where the coefficients \f$ \alpha \f$ and \f$ u_0 \f$ are obtained from the + * \c config, and \f$ u \f$ is the unknown to which the boundary condition is + * applied. + * + * The value \f$ \alpha \cdot [ u_0 - u(x) ] \f$ is a flux. It replaces the + * integrand in the boundary integral for the variable \f$ u \f$. + */ +std::unique_ptr<UniformRobinBoundaryCondition> +createUniformRobinBoundaryCondition( + BaseLib::ConfigTree const& config, + std::vector<MeshLib::Element*>&& elements, + NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id, + int const component_id, unsigned const integration_order, + unsigned const global_dim); + +} // ProcessLib + +#endif // PROCESSLIB_UNIFORMROBINBOUNDARYCONDITION_H diff --git a/ProcessLib/BoundaryCondition/UniformRobinBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/UniformRobinBoundaryConditionLocalAssembler.h new file mode 100644 index 0000000000000000000000000000000000000000..d5d46e907201ed903c30fc3e624f269a52ce8bfc --- /dev/null +++ b/ProcessLib/BoundaryCondition/UniformRobinBoundaryConditionLocalAssembler.h @@ -0,0 +1,86 @@ +/** + * \copyright + * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#ifndef PROCESSLIB_UNIFORMROBINBOUNDARYCONDITIONLOCALASSEMBLER_H +#define PROCESSLIB_UNIFORMROBINBOUNDARYCONDITIONLOCALASSEMBLER_H + +#include "NumLib/DOF/DOFTableUtil.h" + +namespace ProcessLib +{ +struct UniformRobinBoundaryConditionData final { + double const alpha; + double const u_0; +}; + +template <typename ShapeFunction, typename IntegrationMethod, + unsigned GlobalDim> +class UniformRobinBoundaryConditionLocalAssembler final + : public GenericNaturalBoundaryConditionLocalAssembler< + ShapeFunction, IntegrationMethod, GlobalDim> +{ + using Base = GenericNaturalBoundaryConditionLocalAssembler< + ShapeFunction, IntegrationMethod, GlobalDim>; + +public: + UniformRobinBoundaryConditionLocalAssembler( + MeshLib::Element const& e, std::size_t const local_matrix_size, + unsigned const integration_order, + UniformRobinBoundaryConditionData const& data) + : Base(e, integration_order), + _data(data), + _local_K(local_matrix_size, local_matrix_size), + _local_rhs(local_matrix_size) + { + } + + // TODO also implement derivative for Jacobian in Newton scheme. + void assemble(std::size_t const id, + NumLib::LocalToGlobalIndexMap const& dof_table_boundary, + double const t, const GlobalVector& /*x*/, GlobalMatrix& K, + GlobalVector& b) override + { + (void)t; // TODO time-dependent Robin BCs + + _local_K.setZero(); + _local_rhs.setZero(); + + IntegrationMethod integration_method(Base::_integration_order); + std::size_t const n_integration_points = + integration_method.getNumberOfPoints(); + + for (std::size_t ip = 0; ip < n_integration_points; ++ip) { + auto const& sm = Base::_shape_matrices[ip]; + auto const& wp = integration_method.getWeightedPoint(ip); + + // flux = alpha * ( u_0 - u ) + // adding a alpha term to the diagonal of the stiffness matrix + // and a alpha * u_0 term to the rhs vector + _local_K.diagonal().noalias() += + sm.N * _data.alpha * sm.detJ * wp.getWeight(); + _local_rhs.noalias() += + sm.N * _data.alpha * _data.u_0 * sm.detJ * wp.getWeight(); + } + + auto const indices = NumLib::getIndices(id, dof_table_boundary); + K.add(NumLib::LocalToGlobalIndexMap::RowColumnIndices(indices, indices), + _local_K); + b.add(indices, _local_rhs); + } + +private: + UniformRobinBoundaryConditionData const& _data; + + typename Base::NodalMatrixType _local_K; + typename Base::NodalVectorType _local_rhs; +}; + +} // ProcessLib + +#endif // PROCESSLIB_UNIFORMROBINBOUNDARYCONDITIONLOCALASSEMBLER_H diff --git a/Tests/Data b/Tests/Data index 587ff3a5cb2ab33f2923f801242222fa37e46b6e..2aafa43239d3976481dbf0b34574fdd83fa38a12 160000 --- a/Tests/Data +++ b/Tests/Data @@ -1 +1 @@ -Subproject commit 587ff3a5cb2ab33f2923f801242222fa37e46b6e +Subproject commit 2aafa43239d3976481dbf0b34574fdd83fa38a12