From a42275a7e8ff38eb3f0272cb957429fb86286e87 Mon Sep 17 00:00:00 2001 From: Christoph Lehmann Date: Wed, 3 Aug 2016 15:48:40 +0200 Subject: [PATCH 1/8] [PL] registered Robin BC creator --- .../BoundaryCondition/BoundaryCondition.cpp | 39 +++++++++++++------ 1 file changed, 28 insertions(+), 11 deletions(-) diff --git a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp index 1678292f7de..eb182112a80 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" + +std::vector getClonedElements( + MeshGeoToolsLib::BoundaryElementsSearcher& boundary_element_searcher, + GeoLib::GeoObject const& geometry) +{ + std::vector 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 createBoundaryCondition( } else if (type == "UniformNeumann") { - std::vector 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 { From a9ebd7f0c7d9727b067380057a194a11546ca1a2 Mon Sep 17 00:00:00 2001 From: Christoph Lehmann Date: Wed, 3 Aug 2016 15:53:01 +0200 Subject: [PATCH 2/8] [PL] added Robin BC --- .../UniformRobinBoundaryCondition.cpp | 40 +++++++++ .../UniformRobinBoundaryCondition.h | 43 ++++++++++ ...formRobinBoundaryConditionLocalAssembler.h | 86 +++++++++++++++++++ 3 files changed, 169 insertions(+) create mode 100644 ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.cpp create mode 100644 ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.h create mode 100644 ProcessLib/BoundaryCondition/UniformRobinBoundaryConditionLocalAssembler.h diff --git a/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.cpp new file mode 100644 index 00000000000..402bf7633b6 --- /dev/null +++ b/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.cpp @@ -0,0 +1,40 @@ +/** + * \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" +#include "BoundaryConditionConfig.h" +#include "MeshGeoToolsLib/BoundaryElementsSearcher.h" + +namespace ProcessLib +{ +std::unique_ptr +createUniformRobinBoundaryCondition( + BoundaryConditionConfig const& config, + MeshGeoToolsLib::BoundaryElementsSearcher& searcher, + NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id, + unsigned const integration_order, const unsigned global_dim) +{ + DBUG("Constructing RobinBcConfig from config."); + //! \ogs_file_param{boundary_condition__type} + config.config.checkConfigParameter("type", "UniformRobin"); + + //! \ogs_file_param{boundary_condition__UniformRobin__alpha} + double const alpha = config.config.getConfigParameter("alpha"); + //! \ogs_file_param{boundary_condition__UniformRobin__f_0} + double const f_0 = config.config.getConfigParameter("f_0"); + + auto& elems = searcher.getBoundaryElements(config.geometry); + + return std::unique_ptr( + new UniformRobinBoundaryCondition( + integration_order, dof_table, variable_id, config.component_id, + global_dim, elems, UniformRobinBoundaryConditionData{alpha, f_0})); +} + +} // ProcessLib diff --git a/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.h b/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.h new file mode 100644 index 00000000000..91693b9d786 --- /dev/null +++ b/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.h @@ -0,0 +1,43 @@ +/** + * \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$ \partial f/\partial n = \alpha \cdot [ f(x) - f_0 ] \f$, + * where the coefficients \f$ \alpha \f$ and \f$ f_0 \f$ are obtained from the + * \c config. + */ +std::unique_ptr +createUniformRobinBoundaryCondition( + BoundaryConditionConfig const& config, + MeshGeoToolsLib::BoundaryElementsSearcher& searcher, + NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_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 00000000000..13f25c1721e --- /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 f_0; +}; + +template +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); + + // df/dn = alpha ( f - f_0 ) + // adding a -alpha term to the diagonal of the stiffness matrix + // and a -alpha * f_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.f_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 From 54c12bd7e54e56c84b1fc28a3ca1b5051984b92d Mon Sep 17 00:00:00 2001 From: Christoph Lehmann Date: Wed, 3 Aug 2016 15:54:11 +0200 Subject: [PATCH 3/8] [App] added Test for Robin BC --- Applications/CLI/Tests.cmake | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/Applications/CLI/Tests.cmake b/Applications/CLI/Tests.cmake index 78616519c26..8a912e19c02 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 From f0e066d3407b14d29203128d26bd609c1a72787b Mon Sep 17 00:00:00 2001 From: Christoph Lehmann Date: Tue, 9 Aug 2016 10:31:25 +0200 Subject: [PATCH 4/8] [PL] changed construction of Robin BCs --- .../UniformRobinBoundaryCondition.cpp | 22 +++++++++---------- .../UniformRobinBoundaryCondition.h | 7 +++--- 2 files changed, 14 insertions(+), 15 deletions(-) diff --git a/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.cpp index 402bf7633b6..6f9a0102f7e 100644 --- a/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.cpp +++ b/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.cpp @@ -8,33 +8,31 @@ */ #include "UniformRobinBoundaryCondition.h" -#include "BoundaryConditionConfig.h" -#include "MeshGeoToolsLib/BoundaryElementsSearcher.h" namespace ProcessLib { std::unique_ptr createUniformRobinBoundaryCondition( - BoundaryConditionConfig const& config, - MeshGeoToolsLib::BoundaryElementsSearcher& searcher, + BaseLib::ConfigTree const& config, + std::vector&& elements, NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id, - unsigned const integration_order, const unsigned global_dim) + int const component_id, unsigned const integration_order, + unsigned const global_dim) { DBUG("Constructing RobinBcConfig from config."); //! \ogs_file_param{boundary_condition__type} - config.config.checkConfigParameter("type", "UniformRobin"); + config.checkConfigParameter("type", "UniformRobin"); //! \ogs_file_param{boundary_condition__UniformRobin__alpha} - double const alpha = config.config.getConfigParameter("alpha"); + double const alpha = config.getConfigParameter("alpha"); //! \ogs_file_param{boundary_condition__UniformRobin__f_0} - double const f_0 = config.config.getConfigParameter("f_0"); - - auto& elems = searcher.getBoundaryElements(config.geometry); + double const f_0 = config.getConfigParameter("f_0"); return std::unique_ptr( new UniformRobinBoundaryCondition( - integration_order, dof_table, variable_id, config.component_id, - global_dim, elems, UniformRobinBoundaryConditionData{alpha, f_0})); + integration_order, dof_table, variable_id, component_id, global_dim, + std::move(elements), + UniformRobinBoundaryConditionData{alpha, f_0})); } } // ProcessLib diff --git a/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.h b/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.h index 91693b9d786..91ec8ddff0a 100644 --- a/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.h +++ b/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.h @@ -33,10 +33,11 @@ using UniformRobinBoundaryCondition = GenericNaturalBoundaryCondition< */ std::unique_ptr createUniformRobinBoundaryCondition( - BoundaryConditionConfig const& config, - MeshGeoToolsLib::BoundaryElementsSearcher& searcher, + BaseLib::ConfigTree const& config, + std::vector&& elements, NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id, - unsigned const integration_order, unsigned const global_dim); + int const component_id, unsigned const integration_order, + unsigned const global_dim); } // ProcessLib From 78c45ff6b7a29ae57f14109e50ca67003434b1af Mon Sep 17 00:00:00 2001 From: Christoph Lehmann Date: Tue, 9 Aug 2016 15:48:19 +0200 Subject: [PATCH 5/8] [Data] added Robin BC tests --- Tests/Data | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Tests/Data b/Tests/Data index 587ff3a5cb2..8e8c6c93b91 160000 --- a/Tests/Data +++ b/Tests/Data @@ -1 +1 @@ -Subproject commit 587ff3a5cb2ab33f2923f801242222fa37e46b6e +Subproject commit 8e8c6c93b9198edb3588f1c098b1cce1e143b28b From ad204675f2f892d00c821419d82ab1e22e95e7ba Mon Sep 17 00:00:00 2001 From: Christoph Lehmann Date: Wed, 10 Aug 2016 16:40:38 +0200 Subject: [PATCH 6/8] [PL] f-f_0 -> u_0 - u; fixed docu. --- .../UniformRobinBoundaryCondition.cpp | 6 +++--- .../UniformRobinBoundaryCondition.h | 10 +++++++--- .../UniformRobinBoundaryConditionLocalAssembler.h | 14 +++++++------- 3 files changed, 17 insertions(+), 13 deletions(-) diff --git a/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.cpp index 6f9a0102f7e..95bc0c90ca5 100644 --- a/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.cpp +++ b/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.cpp @@ -25,14 +25,14 @@ createUniformRobinBoundaryCondition( //! \ogs_file_param{boundary_condition__UniformRobin__alpha} double const alpha = config.getConfigParameter("alpha"); - //! \ogs_file_param{boundary_condition__UniformRobin__f_0} - double const f_0 = config.getConfigParameter("f_0"); + //! \ogs_file_param{boundary_condition__UniformRobin__u_0} + double const u_0 = config.getConfigParameter("u_0"); return std::unique_ptr( new UniformRobinBoundaryCondition( integration_order, dof_table, variable_id, component_id, global_dim, std::move(elements), - UniformRobinBoundaryConditionData{alpha, f_0})); + UniformRobinBoundaryConditionData{alpha, u_0})); } } // ProcessLib diff --git a/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.h b/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.h index 91ec8ddff0a..87ef77b8825 100644 --- a/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.h +++ b/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.h @@ -27,9 +27,13 @@ using UniformRobinBoundaryCondition = GenericNaturalBoundaryCondition< /*! Creates a new uniform Robin boundary condition from the given data. * * The Robin boundary condition is given in the form - * \f$ \partial f/\partial n = \alpha \cdot [ f(x) - f_0 ] \f$, - * where the coefficients \f$ \alpha \f$ and \f$ f_0 \f$ are obtained from the - * \c config. + * \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 createUniformRobinBoundaryCondition( diff --git a/ProcessLib/BoundaryCondition/UniformRobinBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/UniformRobinBoundaryConditionLocalAssembler.h index 13f25c1721e..d5d46e90720 100644 --- a/ProcessLib/BoundaryCondition/UniformRobinBoundaryConditionLocalAssembler.h +++ b/ProcessLib/BoundaryCondition/UniformRobinBoundaryConditionLocalAssembler.h @@ -16,7 +16,7 @@ namespace ProcessLib { struct UniformRobinBoundaryConditionData final { double const alpha; - double const f_0; + double const u_0; }; template Date: Wed, 10 Aug 2016 16:41:14 +0200 Subject: [PATCH 7/8] [Data] adapted to changes in Robin BC --- Tests/Data | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Tests/Data b/Tests/Data index 8e8c6c93b91..2aafa43239d 160000 --- a/Tests/Data +++ b/Tests/Data @@ -1 +1 @@ -Subproject commit 8e8c6c93b9198edb3588f1c098b1cce1e143b28b +Subproject commit 2aafa43239d3976481dbf0b34574fdd83fa38a12 From a2836b5c054458fe748e23a2e54ccbf8773821d2 Mon Sep 17 00:00:00 2001 From: Christoph Lehmann Date: Wed, 10 Aug 2016 16:49:33 +0200 Subject: [PATCH 8/8] [PL] added "static" --- ProcessLib/BoundaryCondition/BoundaryCondition.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp index eb182112a80..65ce716963a 100644 --- a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp +++ b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp @@ -15,7 +15,7 @@ #include "UniformNeumannBoundaryCondition.h" #include "UniformRobinBoundaryCondition.h" -std::vector getClonedElements( +static std::vector getClonedElements( MeshGeoToolsLib::BoundaryElementsSearcher& boundary_element_searcher, GeoLib::GeoObject const& geometry) {