diff --git a/Applications/DataHolderLib/BoundaryCondition.cpp b/Applications/DataHolderLib/BoundaryCondition.cpp index 47fd4726e9c7bc708d2d2e16dbc24b6a8db0349d..c5b4abbae49fc9123c6c9276baaa109e54acaa38 100644 --- a/Applications/DataHolderLib/BoundaryCondition.cpp +++ b/Applications/DataHolderLib/BoundaryCondition.cpp @@ -28,18 +28,10 @@ BoundaryCondition::ConditionType BoundaryCondition::convertStringToType( { return ConditionType::DIRICHLET; } - if (str == "NonuniformDirichlet") - { - return ConditionType::NONUNIFORMDIRICHLET; - } else if (str == "Neumann") { return ConditionType::NEUMANN; } - else if (str == "NonuniformNeumann") - { - return ConditionType::NONUNIFORMNEUMANN; - } else if (str == "Robin") { return ConditionType::ROBIN; @@ -54,18 +46,10 @@ std::string BoundaryCondition::convertTypeToString(ConditionType type) { return "Dirichlet"; } - if (type == ConditionType::NONUNIFORMDIRICHLET) - { - return "NonuniformDirichlet"; - } else if (type == ConditionType::NEUMANN) { return "Neumann"; } - else if (type == ConditionType::NONUNIFORMNEUMANN) - { - return "NonuniformNeumann"; - } else if (type == ConditionType::ROBIN) { return "Robin"; diff --git a/Applications/DataHolderLib/BoundaryCondition.h b/Applications/DataHolderLib/BoundaryCondition.h index 7ce7f1e061b0e1e9c1ae509ffae4c2030bdeec3f..077ff996fa6a8a230ebc1d74c1f8c57bca94de73 100644 --- a/Applications/DataHolderLib/BoundaryCondition.h +++ b/Applications/DataHolderLib/BoundaryCondition.h @@ -21,9 +21,7 @@ public: { NONE, DIRICHLET, - NONUNIFORMDIRICHLET, NEUMANN, - NONUNIFORMNEUMANN, ROBIN }; diff --git a/Applications/FileIO/XmlIO/OpenGeoSysProject.xsd b/Applications/FileIO/XmlIO/OpenGeoSysProject.xsd index ac0de797669e20f90c15ec11c9040fe1694564b2..632db0d229867441db0e3d1409020e6307d26c89 100644 --- a/Applications/FileIO/XmlIO/OpenGeoSysProject.xsd +++ b/Applications/FileIO/XmlIO/OpenGeoSysProject.xsd @@ -27,6 +27,7 @@ <xs:sequence> <xs:element ref="name" minOccurs="1" maxOccurs="1" /> <xs:element name="type" type="xs:string" minOccurs="1" maxOccurs="1" /> + <xs:element name="mesh" type="xs:string" minOccurs="0" maxOccurs="1" /> <xs:element name="value" minOccurs="0" maxOccurs="1" /> <xs:element name="field_name" type="xs:string" minOccurs="0" maxOccurs="1" /> </xs:sequence> @@ -38,7 +39,6 @@ <xs:element name="geometry" type="xs:string" minOccurs="0" /> <xs:element name="type" type="xs:string" /> <xs:element name="mesh" type="xs:string" minOccurs="0" /> - <xs:element name="field_name" type="xs:string" minOccurs="0" /> <xs:element name="parameter" type="xs:string" minOccurs="0" /> </xs:sequence> </xs:complexType> diff --git a/Documentation/ProjectFile/prj/parameters/parameter/t_mesh.md b/Documentation/ProjectFile/prj/parameters/parameter/t_mesh.md new file mode 100644 index 0000000000000000000000000000000000000000..65e9f2aca8d34f125ce551aba0fa2cc648a70854 --- /dev/null +++ b/Documentation/ProjectFile/prj/parameters/parameter/t_mesh.md @@ -0,0 +1,6 @@ +The domain of definition of the parameter. + +The parameter's domain of definition is implicitly the main mesh. +It needs to be specified explicitly for parameters defined on boundary +or subdomain meshes. A notable exception is the Constant parameter, +which has an arbitrary domain of definition. diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformDirichlet/c_NonuniformDirichlet.md b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformDirichlet/c_NonuniformDirichlet.md deleted file mode 100644 index a4df55e6ec1007021690f91c9933f42783d21ff5..0000000000000000000000000000000000000000 --- a/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformDirichlet/c_NonuniformDirichlet.md +++ /dev/null @@ -1,2 +0,0 @@ -Declares a Dirichlet boundary condition that is non-uniform in space. -The values are read from the given boundary mesh. diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformDirichlet/t_field_name.md b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformDirichlet/t_field_name.md deleted file mode 100644 index fc5aead79e2a9cd6e2b973505bc4c317288903e4..0000000000000000000000000000000000000000 --- a/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformDirichlet/t_field_name.md +++ /dev/null @@ -1,3 +0,0 @@ -This specifies a nodal field defined on the surface mesh (see the \c mesh -parameter of the nonuniform Dirichlet BC) and the values used for the boundary -condition. diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformNeumann/c_NonuniformNeumann.md b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformNeumann/c_NonuniformNeumann.md deleted file mode 100644 index cf46bf6546ca6b7816fb337f8920904f8173952c..0000000000000000000000000000000000000000 --- a/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformNeumann/c_NonuniformNeumann.md +++ /dev/null @@ -1 +0,0 @@ -Declares a Neumann boundary condition that is non-uniform in space. diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformNeumann/t_field_name.md b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformNeumann/t_field_name.md deleted file mode 100644 index 7f36683addd9a72514cca3010606b6cd24bab004..0000000000000000000000000000000000000000 --- a/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformNeumann/t_field_name.md +++ /dev/null @@ -1,5 +0,0 @@ -This specifies a nodal field defined on the surface mesh (see the \c mesh -parameter of the nonuniform Neumann BC). - -From that nodal field the flux values, which are applied via this Neumann BC, -are read. diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformVariableDependentNeumann/c_NonuniformVariableDependentNeumann.md b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/VariableDependentNeumann/c_NonuniformVariableDependentNeumann.md similarity index 100% rename from Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformVariableDependentNeumann/c_NonuniformVariableDependentNeumann.md rename to Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/VariableDependentNeumann/c_NonuniformVariableDependentNeumann.md diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformVariableDependentNeumann/t_coefficient_current_variable_name.md b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/VariableDependentNeumann/t_coefficient_current_variable_name.md similarity index 100% rename from Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformVariableDependentNeumann/t_coefficient_current_variable_name.md rename to Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/VariableDependentNeumann/t_coefficient_current_variable_name.md diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformVariableDependentNeumann/t_coefficient_mixed_variables_name.md b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/VariableDependentNeumann/t_coefficient_mixed_variables_name.md similarity index 100% rename from Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformVariableDependentNeumann/t_coefficient_mixed_variables_name.md rename to Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/VariableDependentNeumann/t_coefficient_mixed_variables_name.md diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformVariableDependentNeumann/t_coefficient_other_variable_name.md b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/VariableDependentNeumann/t_coefficient_other_variable_name.md similarity index 100% rename from Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformVariableDependentNeumann/t_coefficient_other_variable_name.md rename to Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/VariableDependentNeumann/t_coefficient_other_variable_name.md diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformVariableDependentNeumann/t_constant_name.md b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/VariableDependentNeumann/t_constant_name.md similarity index 100% rename from Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformVariableDependentNeumann/t_constant_name.md rename to Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/VariableDependentNeumann/t_constant_name.md diff --git a/MeshLib/Mesh.h b/MeshLib/Mesh.h index 44d0059239ed2f5cec76dcc62a56f4278830e7b2..5cf9460d5dd9f247a1d1887f19f708f105c8cfa1 100644 --- a/MeshLib/Mesh.h +++ b/MeshLib/Mesh.h @@ -188,6 +188,16 @@ protected: }; /* class */ +/// Meshes are equal if their id's are equal. +inline bool operator==(Mesh const& a, Mesh const& b) +{ + return a.getID() == b.getID(); +} + +inline bool operator!=(Mesh const& a, Mesh const& b) +{ + return !(a == b); +} /// Scales the mesh property with name \c property_name by given \c factor. /// \note The property must be a "double" property. diff --git a/ProcessLib/BoundaryCondition/CreateBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/CreateBoundaryCondition.cpp index e3518dd4c0b5b8edc069192d506f336c9f93b088..4efc1805feff70ab82ca80a24de38952e3bf9bdd 100644 --- a/ProcessLib/BoundaryCondition/CreateBoundaryCondition.cpp +++ b/ProcessLib/BoundaryCondition/CreateBoundaryCondition.cpp @@ -15,12 +15,10 @@ #include "DirichletBoundaryCondition.h" #include "DirichletBoundaryConditionWithinTimeInterval.h" #include "NeumannBoundaryCondition.h" -#include "NonuniformDirichletBoundaryCondition.h" -#include "NonuniformNeumannBoundaryCondition.h" -#include "NonuniformVariableDependentNeumannBoundaryCondition.h" #include "NormalTractionBoundaryCondition.h" #include "PhaseFieldIrreversibleDamageOracleBoundaryCondition.h" #include "RobinBoundaryCondition.h" +#include "VariableDependentNeumannBoundaryCondition.h" #include "BaseLib/TimeInterval.h" @@ -78,27 +76,12 @@ std::unique_ptr<BoundaryCondition> createBoundaryCondition( *config.component_id, integration_order, shapefunction_order, bulk_mesh.getDimension(), parameters); } - if (type == "NonuniformDirichlet") + if (type == "VariableDependentNeumann") { - return ProcessLib::createNonuniformDirichletBoundaryCondition( - config.config, config.boundary_mesh, dof_table, variable_id, - *config.component_id); - } - if (type == "NonuniformNeumann") - { - return ProcessLib::createNonuniformNeumannBoundaryCondition( + return ProcessLib::createVariableDependentNeumannBoundaryCondition( config.config, config.boundary_mesh, dof_table, variable_id, *config.component_id, integration_order, shapefunction_order, - bulk_mesh); - } - - if (type == "NonuniformVariableDependentNeumann") - { - return ProcessLib:: - createNonuniformVariableDependentNeumannBoundaryCondition( - config.config, config.boundary_mesh, dof_table, variable_id, - *config.component_id, integration_order, shapefunction_order, - bulk_mesh); + bulk_mesh.getDimension(), parameters); } if (type == "Python") diff --git a/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.cpp index fc757a7b5c0db70d132d9622b4dca63c233d8539..c714a531dec83f6ec4c75e076bfd18b546aa4d0d 100644 --- a/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.cpp +++ b/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.cpp @@ -67,7 +67,16 @@ std::unique_ptr<DirichletBoundaryCondition> createDirichletBoundaryCondition( auto const param_name = config.getConfigParameter<std::string>("parameter"); DBUG("Using parameter %s", param_name.c_str()); - auto& param = findParameter<double>(param_name, parameters, 1); + auto& parameter = findParameter<double>(param_name, parameters, 1); + + if (parameter.mesh() && *parameter.mesh() != bc_mesh) + { + OGS_FATAL( + "The boundary condition is defined on a different domain than the " + "parameter: boundary condition is defined on mesh '%s', parameter " + "is defined on mesh '%s'.", + bc_mesh.getName().c_str(), parameter.mesh()->getName().c_str()); + } // In case of partitioned mesh the boundary could be empty, i.e. there is no // boundary condition. @@ -84,7 +93,7 @@ std::unique_ptr<DirichletBoundaryCondition> createDirichletBoundaryCondition( #endif // USE_PETSC return std::make_unique<DirichletBoundaryCondition>( - param, bc_mesh, dof_table_bulk, variable_id, component_id); + parameter, bc_mesh, dof_table_bulk, variable_id, component_id); } } // namespace ProcessLib diff --git a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryConditionLocalAssembler.h index e0c086be69b49f61b365eddc79b4a278fad46705..12fbf69dbd7ffd45b2fc72bc10e2bcbef88c0a62 100644 --- a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryConditionLocalAssembler.h +++ b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryConditionLocalAssembler.h @@ -38,24 +38,58 @@ protected: using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType; using NodalVectorType = typename ShapeMatricesType::NodalVectorType; + struct NAndWeight + { + NAndWeight(typename ShapeMatricesType::ShapeMatrices::ShapeType&& N_, + double const weight_) + : N(std::move(N_)), weight(weight_) + { + } + typename ShapeMatricesType::ShapeMatrices::ShapeType const N; + double const weight; + }; + +private: + static std::vector<NAndWeight, Eigen::aligned_allocator<NAndWeight>> + initNsAndWeights(MeshLib::Element const& e, bool is_axially_symmetric, + unsigned const integration_order) + { + IntegrationMethod const integration_method(integration_order); + std::vector<NAndWeight, Eigen::aligned_allocator<NAndWeight>> + ns_and_weights; + ns_and_weights.reserve(integration_method.getNumberOfPoints()); + + auto sms = initShapeMatrices<ShapeFunction, ShapeMatricesType, + IntegrationMethod, GlobalDim>( + e, is_axially_symmetric, integration_method); + for (unsigned ip = 0; ip < sms.size(); ++ip) + { + auto& sm = sms[ip]; + double const w = + sm.detJ * sm.integralMeasure * + integration_method.getWeightedPoint(ip).getWeight(); + + ns_and_weights.emplace_back(std::move(sm.N), w); + } + + return ns_and_weights; + } + public: GenericNaturalBoundaryConditionLocalAssembler( MeshLib::Element const& e, bool is_axially_symmetric, unsigned const integration_order) : _integration_method(integration_order), - _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType, - IntegrationMethod, GlobalDim>( - e, is_axially_symmetric, _integration_method)), + _ns_and_weights( + initNsAndWeights(e, is_axially_symmetric, integration_order)), _element(e) { } protected: IntegrationMethod const _integration_method; - std::vector<typename ShapeMatricesType::ShapeMatrices, - Eigen::aligned_allocator< - typename ShapeMatricesType::ShapeMatrices>> const - _shape_matrices; + std::vector<NAndWeight, Eigen::aligned_allocator<NAndWeight>> const + _ns_and_weights; MeshLib::Element const& _element; }; diff --git a/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition-impl.h b/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition-impl.h deleted file mode 100644 index 39ac47c9dad5fccd5072a010a371215941ef0838..0000000000000000000000000000000000000000 --- a/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition-impl.h +++ /dev/null @@ -1,87 +0,0 @@ -/** - * \copyright - * Copyright (c) 2012-2019, 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 "GenericNonuniformNaturalBoundaryConditionLocalAssembler.h" -#include "ProcessLib/Utils/CreateLocalAssemblers.h" - -namespace ProcessLib -{ -template <typename BoundaryConditionData, - template <typename, typename, unsigned> - class LocalAssemblerImplementation> -template <typename Data> -GenericNonuniformNaturalBoundaryCondition<BoundaryConditionData, - LocalAssemblerImplementation>:: - GenericNonuniformNaturalBoundaryCondition( - unsigned const integration_order, unsigned const shapefunction_order, - NumLib::LocalToGlobalIndexMap const& dof_table_bulk, - int const variable_id, int const component_id, - unsigned const global_dim, MeshLib::Mesh const& bc_mesh, Data&& data) - : _data(std::forward<Data>(data)), _bc_mesh(bc_mesh) -{ - static_assert(std::is_same<typename std::decay<BoundaryConditionData>::type, - typename std::decay<Data>::type>::value, - "Type mismatch between declared and passed BC data."); - - // check basic data consistency - if (variable_id >= - static_cast<int>(dof_table_bulk.getNumberOfVariables()) || - component_id >= - dof_table_bulk.getNumberOfVariableComponents(variable_id)) - { - OGS_FATAL( - "Variable id or component id too high. Actual values: (%d, %d), " - "maximum values: (%d, %d).", - variable_id, component_id, dof_table_bulk.getNumberOfVariables(), - dof_table_bulk.getNumberOfVariableComponents(variable_id)); - } - - if (_bc_mesh.getDimension() + 1 != global_dim) - { - OGS_FATAL( - "The dimension of the given boundary mesh (%d) is not by one lower " - "than the bulk dimension (%d).", - _bc_mesh.getDimension(), global_dim); - } - - std::vector<MeshLib::Node*> const& bc_nodes = _bc_mesh.getNodes(); - DBUG("Found %d nodes for Natural BCs for the variable %d and component %d", - bc_nodes.size(), variable_id, component_id); - - MeshLib::MeshSubset bc_mesh_subset(_bc_mesh, bc_nodes); - - // Create local DOF table from the BC mesh subset for the given variable and - // component id. - _dof_table_boundary.reset(dof_table_bulk.deriveBoundaryConstrainedMap( - variable_id, {component_id}, std::move(bc_mesh_subset))); - - createLocalAssemblers<LocalAssemblerImplementation>( - global_dim, _bc_mesh.getElements(), *_dof_table_boundary, - shapefunction_order, _local_assemblers, _bc_mesh.isAxiallySymmetric(), - integration_order, _data); -} - -template <typename BoundaryConditionData, - template <typename, typename, unsigned> - class LocalAssemblerImplementation> -void GenericNonuniformNaturalBoundaryCondition< - BoundaryConditionData, - LocalAssemblerImplementation>::applyNaturalBC(const double t, - const GlobalVector& x, - GlobalMatrix& K, - GlobalVector& b, - GlobalMatrix* Jac) -{ - GlobalExecutor::executeMemberOnDereferenced( - &GenericNonuniformNaturalBoundaryConditionLocalAssemblerInterface:: - assemble, - _local_assemblers, *_dof_table_boundary, t, x, K, b, Jac); -} - -} // namespace ProcessLib diff --git a/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition.h b/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition.h deleted file mode 100644 index f308222367490039c62ee42620e1673078bc06fc..0000000000000000000000000000000000000000 --- a/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition.h +++ /dev/null @@ -1,59 +0,0 @@ -/** - * \copyright - * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org) - * Distributed under a Modified BSD License. - * See accompanying file LICENSE.txt or - * http://www.opengeosys.org/project/license - * - */ - -#pragma once - -#include "BoundaryCondition.h" -#include "MeshLib/MeshSubset.h" - -namespace ProcessLib -{ -class GenericNonuniformNaturalBoundaryConditionLocalAssemblerInterface; - -template <typename BoundaryConditionData, - template <typename, typename, unsigned> - class LocalAssemblerImplementation> -class GenericNonuniformNaturalBoundaryCondition final : public BoundaryCondition -{ -public: - /// Create a boundary condition process from given config, - /// DOF-table, and a mesh subset for a given variable and its component. - /// A local DOF-table, a subset of the given one, is constructed. - template <typename Data> - GenericNonuniformNaturalBoundaryCondition( - unsigned const integration_order, unsigned const shapefunction_order, - NumLib::LocalToGlobalIndexMap const& dof_table_bulk, - int const variable_id, int const component_id, - unsigned const global_dim, MeshLib::Mesh const& bc_mesh, Data&& data); - - /// Calls local assemblers which calculate their contributions to the global - /// matrix and the right-hand-side. - void applyNaturalBC(const double t, GlobalVector const& x, GlobalMatrix& K, - GlobalVector& b, GlobalMatrix* Jac) override; - -private: - /// Data used in the assembly of the specific boundary condition. - BoundaryConditionData _data; - - /// A lower-dimensional mesh on which the boundary condition is defined. - MeshLib::Mesh const& _bc_mesh; - - /// Local dof table, a subset of the global one restricted to the - /// participating number of _elements of the boundary condition. - std::unique_ptr<NumLib::LocalToGlobalIndexMap> _dof_table_boundary; - - /// Local assemblers for each element of the boundary mesh. - std::vector<std::unique_ptr< - GenericNonuniformNaturalBoundaryConditionLocalAssemblerInterface>> - _local_assemblers; -}; - -} // namespace ProcessLib - -#include "GenericNonuniformNaturalBoundaryCondition-impl.h" diff --git a/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryConditionLocalAssembler.h deleted file mode 100644 index e6a4add85c5e247c9ce6555bedde39c7044aeb9b..0000000000000000000000000000000000000000 --- a/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryConditionLocalAssembler.h +++ /dev/null @@ -1,96 +0,0 @@ -/** - * \copyright - * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org) - * Distributed under a Modified BSD License. - * See accompanying file LICENSE.txt or - * http://www.opengeosys.org/project/license - * - */ - -#pragma once - -#include "NumLib/DOF/LocalToGlobalIndexMap.h" -#include "NumLib/Fem/ShapeMatrixPolicy.h" -#include "ProcessLib/Utils/InitShapeMatrices.h" - -namespace ProcessLib -{ -class GenericNonuniformNaturalBoundaryConditionLocalAssemblerInterface -{ -public: - virtual ~GenericNonuniformNaturalBoundaryConditionLocalAssemblerInterface() = - default; - - virtual void assemble( - std::size_t const id, - NumLib::LocalToGlobalIndexMap const& dof_table_boundary, double const t, - const GlobalVector& x, GlobalMatrix& K, GlobalVector& b, - GlobalMatrix* Jac) = 0; -}; - -template <typename ShapeFunction, typename IntegrationMethod, - unsigned GlobalDim> -class GenericNonuniformNaturalBoundaryConditionLocalAssembler - : public GenericNonuniformNaturalBoundaryConditionLocalAssemblerInterface -{ -protected: - using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>; - using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType; - using NodalVectorType = typename ShapeMatricesType::NodalVectorType; - - struct NAndWeight - { - NAndWeight(typename ShapeMatricesType::ShapeMatrices::ShapeType&& N_, - double const weight_) - : N(std::move(N_)), weight(weight_) - { - } - typename ShapeMatricesType::ShapeMatrices::ShapeType const N; - double const weight; - }; - -private: - static std::vector<NAndWeight, Eigen::aligned_allocator<NAndWeight>> - initNsAndWeights(MeshLib::Element const& e, bool is_axially_symmetric, - unsigned const integration_order) - { - IntegrationMethod const integration_method(integration_order); - std::vector<NAndWeight, Eigen::aligned_allocator<NAndWeight>> - ns_and_weights; - ns_and_weights.reserve(integration_method.getNumberOfPoints()); - - auto sms = initShapeMatrices<ShapeFunction, ShapeMatricesType, - IntegrationMethod, GlobalDim>( - e, is_axially_symmetric, integration_method); - for (unsigned ip = 0; ip < sms.size(); ++ip) - { - auto& sm = sms[ip]; - double const w = - sm.detJ * sm.integralMeasure * - integration_method.getWeightedPoint(ip).getWeight(); - - ns_and_weights.emplace_back(std::move(sm.N), w); - } - - return ns_and_weights; - } - -public: - GenericNonuniformNaturalBoundaryConditionLocalAssembler( - MeshLib::Element const& e, bool is_axially_symmetric, - unsigned const integration_order) - : _integration_method(integration_order), - _element(e), - _ns_and_weights( - initNsAndWeights(e, is_axially_symmetric, integration_order)) - { - } - -protected: - IntegrationMethod const _integration_method; - MeshLib::Element const& _element; - std::vector<NAndWeight, Eigen::aligned_allocator<NAndWeight>> const - _ns_and_weights; -}; - -} // namespace ProcessLib diff --git a/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h index 529313734041734872b3440f6a7e5bf55632ec94..0f09c01d8a8d7d5be1a4ee471ef20da288defdfd 100644 --- a/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h +++ b/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h @@ -26,6 +26,7 @@ class NeumannBoundaryConditionLocalAssembler final using Base = GenericNaturalBoundaryConditionLocalAssembler< ShapeFunction, IntegrationMethod, GlobalDim>; using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>; + using NodalVectorType = typename Base::NodalVectorType; public: /// The neumann_bc_value factor is directly integrated into the local @@ -53,26 +54,18 @@ public: unsigned const n_integration_points = Base::_integration_method.getNumberOfPoints(); - SpatialPosition pos; - pos.setElementID(id); - - using FemType = - NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>; - FemType fe( - static_cast<const typename ShapeFunction::MeshElement&>(Base::_element)); + // Get element nodes for the interpolation from nodes to integration + // point. + NodalVectorType parameter_node_values = + _neumann_bc_parameter.getNodalValuesOnElement(Base::_element, t); for (unsigned ip = 0; ip < n_integration_points; ip++) { - pos.setIntegrationPoint(ip); - auto const& sm = Base::_shape_matrices[ip]; - auto const& wp = Base::_integration_method.getWeightedPoint(ip); - - MathLib::TemplatePoint<double, 3> coords_ip(fe.interpolateCoordinates(sm.N)); - pos.setCoordinates(coords_ip); + auto const& ip_data = Base::_ns_and_weights[ip]; - _local_rhs.noalias() += sm.N * _neumann_bc_parameter(t, pos)[0] * - sm.detJ * wp.getWeight() * - sm.integralMeasure; + _local_rhs.noalias() += ip_data.N * + parameter_node_values.dot(ip_data.N) * + ip_data.weight; } auto const indices = NumLib::getIndices(id, dof_table_boundary); @@ -81,7 +74,7 @@ public: private: Parameter<double> const& _neumann_bc_parameter; - typename Base::NodalVectorType _local_rhs; + NodalVectorType _local_rhs; public: EIGEN_MAKE_ALIGNED_OPERATOR_NEW; diff --git a/ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.cpp deleted file mode 100644 index 5786e0947c6f2c3840e12746fe58fdfcdb1ec01b..0000000000000000000000000000000000000000 --- a/ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.cpp +++ /dev/null @@ -1,63 +0,0 @@ -/** - * \file - * - * \copyright - * Copyright (c) 2012-2019, 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 "NonuniformDirichletBoundaryCondition.h" - -#include "MeshLib/IO/readMeshFromFile.h" -#include "ProcessLib/Utils/ProcessUtils.h" - -namespace ProcessLib -{ -std::unique_ptr<NonuniformDirichletBoundaryCondition> -createNonuniformDirichletBoundaryCondition( - BaseLib::ConfigTree const& config, MeshLib::Mesh const& boundary_mesh, - NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id, - int const component_id) -{ - DBUG("Constructing NonuniformDirichlet BC from config."); - //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__type} - config.checkConfigParameter("type", "NonuniformDirichlet"); - - // TODO finally use ProcessLib::Parameter here - auto const field_name = - //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__NonuniformDirichlet__field_name} - config.getConfigParameter<std::string>("field_name"); - - auto const* const property = - boundary_mesh.getProperties().getPropertyVector<double>( - field_name, MeshLib::MeshItemType::Node, 1); - - if (!property) - { - OGS_FATAL("A property with name `%s' does not exist in `%s'.", - field_name.c_str(), boundary_mesh.getName().c_str()); - } - - // In case of partitioned mesh the boundary could be empty, i.e. there is no - // boundary condition. -#ifdef USE_PETSC - // This can be extracted to createBoundaryCondition() but then the config - // parameters are not read and will cause an error. - // TODO (naumov): Add a function to ConfigTree for skipping the tags of the - // subtree and move the code up in createBoundaryCondition(). - if (boundary_mesh.getDimension() == 0 && - boundary_mesh.getNumberOfNodes() == 0 && - boundary_mesh.getNumberOfElements() == 0) - { - return nullptr; - } -#endif // USE_PETSC - - return std::make_unique<NonuniformDirichletBoundaryCondition>( - boundary_mesh, *property, dof_table, variable_id, component_id); -} - -} // namespace ProcessLib diff --git a/ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.h b/ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.h deleted file mode 100644 index ba67cb886f3236a3089ce0f421a5140486f857f8..0000000000000000000000000000000000000000 --- a/ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.h +++ /dev/null @@ -1,129 +0,0 @@ -/** - * \file - * - * \copyright - * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org) - * Distributed under a Modified BSD License. - * See accompanying file LICENSE.txt or - * http://www.opengeosys.org/project/license - * - */ - -#pragma once - -#include "BoundaryCondition.h" - -#include "MeshLib/PropertyVector.h" -#include "NumLib/DOF/DOFTableUtil.h" -#include "NumLib/DOF/LocalToGlobalIndexMap.h" -#include "NumLib/IndexValueVector.h" -#include "ProcessLib/Parameter/Parameter.h" -#include "ProcessLib/Utils/ProcessUtils.h" - -namespace ProcessLib -{ -class NonuniformDirichletBoundaryCondition final : public BoundaryCondition -{ -public: - NonuniformDirichletBoundaryCondition( - // int const bulk_mesh_dimension, - MeshLib::Mesh const& boundary_mesh, - MeshLib::PropertyVector<double> const& values, - NumLib::LocalToGlobalIndexMap const& dof_table_bulk, - int const variable_id, - int const component_id) - : _values(values), - _boundary_mesh(boundary_mesh), - _variable_id(variable_id), - _component_id(component_id) - { - if (_variable_id >= - static_cast<int>(dof_table_bulk.getNumberOfVariables()) || - _component_id >= - dof_table_bulk.getNumberOfVariableComponents(_variable_id)) - { - OGS_FATAL( - "Variable id or component id too high. Actual values: (%d, " - "%d), " - "maximum values: (%d, %d).", - _variable_id, _component_id, - dof_table_bulk.getNumberOfVariables(), - dof_table_bulk.getNumberOfVariableComponents(_variable_id)); - } - - if (!_boundary_mesh.getProperties().existsPropertyVector<std::size_t>( - "bulk_node_ids")) - { - OGS_FATAL( - "The required bulk node ids map does not exist in the boundary " - "mesh '%s'.", - _boundary_mesh.getName().c_str()); - } - - std::vector<MeshLib::Node*> const& bc_nodes = _boundary_mesh.getNodes(); - DBUG( - "Found %d nodes for Natural BCs for the variable %d and component " - "%d", - bc_nodes.size(), variable_id, component_id); - - MeshLib::MeshSubset boundary_mesh_subset(_boundary_mesh, bc_nodes); - - // Create local DOF table from the BC mesh subset for the given variable - // and component id. - _dof_table_boundary.reset(dof_table_bulk.deriveBoundaryConstrainedMap( - variable_id, {component_id}, std::move(boundary_mesh_subset))); - } - - void getEssentialBCValues( - const double /*t*/, GlobalVector const& /*x*/, - NumLib::IndexValueVector<GlobalIndexType>& bc_values) const override - { - bc_values.ids.clear(); - bc_values.values.clear(); - - // Convert mesh node ids to global index for the given component. - bc_values.ids.reserve(_values.size()); - bc_values.values.reserve(_values.size()); - - // Map boundary dof indices to bulk dof indices and the corresponding - // values. - for (auto const* const node : _boundary_mesh.getNodes()) - { - auto const node_id = node->getID(); - auto const global_index = _dof_table_boundary->getGlobalIndex( - {_boundary_mesh.getID(), MeshLib::MeshItemType::Node, node_id}, - _variable_id, _component_id); - if (global_index == NumLib::MeshComponentMap::nop) - { - continue; - } - // For the DDC approach (e.g. with PETSc option), the negative index - // of global_index means that the entry by that index is a ghost - // one, which should be dropped. Especially for PETSc routines - // MatZeroRows and MatZeroRowsColumns, which are called to apply the - // Dirichlet BC, the negative index is not accepted like other - // matrix or vector PETSc routines. Therefore, the following - // if-condition is applied. - if (global_index >= 0) - { - bc_values.ids.emplace_back(global_index); - bc_values.values.push_back(_values[node_id]); - } - } - } - -private: - MeshLib::PropertyVector<double> const& _values; - MeshLib::Mesh const& _boundary_mesh; - std::unique_ptr<NumLib::LocalToGlobalIndexMap const> _dof_table_boundary; - int const _variable_id; - int const _component_id; -}; - -std::unique_ptr<NonuniformDirichletBoundaryCondition> -createNonuniformDirichletBoundaryCondition( - BaseLib::ConfigTree const& config, MeshLib::Mesh const& boundary_mesh, - NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id, - int const component_id); - -} // namespace ProcessLib diff --git a/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.cpp deleted file mode 100644 index bbe963d3cae71e39349c59936f8ab09d2a3d15c0..0000000000000000000000000000000000000000 --- a/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.cpp +++ /dev/null @@ -1,92 +0,0 @@ -/** - * \copyright - * Copyright (c) 2012-2019, 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 "NonuniformNeumannBoundaryCondition.h" - -#include "MeshLib/IO/readMeshFromFile.h" -#include "ProcessLib/Utils/ProcessUtils.h" - -namespace ProcessLib -{ -std::unique_ptr<NonuniformNeumannBoundaryCondition> -createNonuniformNeumannBoundaryCondition( - BaseLib::ConfigTree const& config, MeshLib::Mesh const& boundary_mesh, - NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id, - int const component_id, unsigned const integration_order, - unsigned const shapefunction_order, MeshLib::Mesh const& bulk_mesh) -{ - DBUG("Constructing NonuniformNeumann BC from config."); - //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__type} - config.checkConfigParameter("type", "NonuniformNeumann"); - - // TODO finally use ProcessLib::Parameter here - auto const field_name = - //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__NonuniformNeumann__field_name} - config.getConfigParameter<std::string>("field_name"); - - auto const* const property = - boundary_mesh.getProperties().getPropertyVector<double>(field_name); - - if (!property) - { - OGS_FATAL("A property with name `%s' does not exist in `%s'.", - field_name.c_str(), boundary_mesh.getName().c_str()); - } - - if (property->getMeshItemType() != MeshLib::MeshItemType::Node) - { - OGS_FATAL( - "Only nodal fields are supported for non-uniform fields. Field " - "`%s' is not nodal.", - field_name.c_str()); - } - - if (property->getNumberOfComponents() != 1) - { - OGS_FATAL("`%s' is not a one-component field.", field_name.c_str()); - } - - std::string const mapping_to_bulk_nodes_property = "bulk_node_ids"; - auto const* const mapping_to_bulk_nodes = - boundary_mesh.getProperties().getPropertyVector<std::size_t>( - mapping_to_bulk_nodes_property); - - if (mapping_to_bulk_nodes == nullptr || - mapping_to_bulk_nodes->getMeshItemType() != - MeshLib::MeshItemType::Node || - mapping_to_bulk_nodes->getNumberOfComponents() != 1) - { - OGS_FATAL("Field `%s' is not set up properly.", - mapping_to_bulk_nodes_property.c_str()); - } - - // In case of partitioned mesh the boundary could be empty, i.e. there is no - // boundary condition. -#ifdef USE_PETSC - // This can be extracted to createBoundaryCondition() but then the config - // parameters are not read and will cause an error. - // TODO (naumov): Add a function to ConfigTree for skipping the tags of the - // subtree and move the code up in createBoundaryCondition(). - if (boundary_mesh.getDimension() == 0 && - boundary_mesh.getNumberOfNodes() == 0 && - boundary_mesh.getNumberOfElements() == 0) - { - return nullptr; - } -#endif // USE_PETSC - - return std::make_unique<NonuniformNeumannBoundaryCondition>( - integration_order, shapefunction_order, dof_table, variable_id, - component_id, bulk_mesh.getDimension(), boundary_mesh, - NonuniformNeumannBoundaryConditionData{ - *property, bulk_mesh.getID(), *mapping_to_bulk_nodes, dof_table, - variable_id, component_id}); -} - -} // namespace ProcessLib diff --git a/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryConditionLocalAssembler.h deleted file mode 100644 index 88ea7ad1686d4fd672dc9e1e0d8eb644d20e5b7b..0000000000000000000000000000000000000000 --- a/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryConditionLocalAssembler.h +++ /dev/null @@ -1,95 +0,0 @@ -/** - * \copyright - * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org) - * Distributed under a Modified BSD License. - * See accompanying file LICENSE.txt or - * http://www.opengeosys.org/project/license - * - */ - -#pragma once - -#include "MeshLib/PropertyVector.h" -#include "NumLib/DOF/DOFTableUtil.h" -#include "NumLib/Function/Interpolation.h" -#include "ProcessLib/Parameter/MeshNodeParameter.h" - -#include "GenericNonuniformNaturalBoundaryConditionLocalAssembler.h" - -namespace ProcessLib -{ -struct NonuniformNeumannBoundaryConditionData -{ - /* TODO use Parameter */ - MeshLib::PropertyVector<double> const& values; - - // Used for mapping boundary nodes to bulk nodes. - std::size_t bulk_mesh_id; - MeshLib::PropertyVector<std::size_t> const& mapping_to_bulk_nodes; - NumLib::LocalToGlobalIndexMap const& dof_table_bulk; - int const variable_id_bulk; - int const component_id_bulk; -}; - -template <typename ShapeFunction, typename IntegrationMethod, - unsigned GlobalDim> -class NonuniformNeumannBoundaryConditionLocalAssembler final - : public GenericNonuniformNaturalBoundaryConditionLocalAssembler< - ShapeFunction, IntegrationMethod, GlobalDim> -{ - using Base = GenericNonuniformNaturalBoundaryConditionLocalAssembler< - ShapeFunction, IntegrationMethod, GlobalDim>; - using NodalVectorType = typename Base::NodalVectorType; - -public: - /// The neumann_bc_value factor is directly integrated into the local - /// element matrix. - NonuniformNeumannBoundaryConditionLocalAssembler( - MeshLib::Element const& e, - std::size_t const local_matrix_size, - bool const is_axially_symmetric, - unsigned const integration_order, - NonuniformNeumannBoundaryConditionData const& data) - : Base(e, is_axially_symmetric, integration_order), - _data(data), - _local_rhs(local_matrix_size) - { - } - - void assemble(std::size_t const id, - NumLib::LocalToGlobalIndexMap const& dof_table_boundary, - double const t, const GlobalVector& /*x*/, - GlobalMatrix& /*K*/, GlobalVector& b, - GlobalMatrix* /*Jac*/) override - { - _local_rhs.setZero(); - - MeshNodeParameter<double> neumann_values{"NeumannValues", _data.values}; - // Get element nodes for the interpolation from nodes to - // integration point. - NodalVectorType parameter_node_values = - neumann_values.getNodalValuesOnElement(Base::_element, t); - - unsigned const n_integration_points = - Base::_integration_method.getNumberOfPoints(); - for (unsigned ip = 0; ip < n_integration_points; ip++) - { - auto const& n_and_weight = Base::_ns_and_weights[ip]; - auto const& N = n_and_weight.N; - auto const& w = n_and_weight.weight; - _local_rhs.noalias() += N * parameter_node_values.dot(N) * w; - } - - auto const indices = NumLib::getIndices(id, dof_table_boundary); - b.add(indices, _local_rhs); - } - -private: - NonuniformNeumannBoundaryConditionData const& _data; - NodalVectorType _local_rhs; - -public: - EIGEN_MAKE_ALIGNED_OPERATOR_NEW; -}; - -} // namespace ProcessLib diff --git a/ProcessLib/BoundaryCondition/NonuniformVariableDependentNeumannBoundaryCondition.h b/ProcessLib/BoundaryCondition/NonuniformVariableDependentNeumannBoundaryCondition.h deleted file mode 100644 index fe7483df0f1fa6ce9ad8540ae3d51df46b5ed9ab..0000000000000000000000000000000000000000 --- a/ProcessLib/BoundaryCondition/NonuniformVariableDependentNeumannBoundaryCondition.h +++ /dev/null @@ -1,30 +0,0 @@ -/** - * \copyright - * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org) - * Distributed under a Modified BSD License. - * See accompanying file LICENSE.txt or - * http://www.opengeosys.org/project/license - * - */ - -#pragma once - -#include "GenericNonuniformNaturalBoundaryCondition.h" -#include "MeshLib/PropertyVector.h" -#include "NonuniformVariableDependentNeumannBoundaryConditionLocalAssembler.h" - -namespace ProcessLib -{ -using NonuniformVariableDependentNeumannBoundaryCondition = - GenericNonuniformNaturalBoundaryCondition< - NonuniformVariableDependentNeumannBoundaryConditionData, - NonuniformVariableDependentNeumannBoundaryConditionLocalAssembler>; - -std::unique_ptr<NonuniformVariableDependentNeumannBoundaryCondition> -createNonuniformVariableDependentNeumannBoundaryCondition( - BaseLib::ConfigTree const& config, MeshLib::Mesh const& boundary_mesh, - NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id, - int const component_id, unsigned const integration_order, - unsigned const shapefunction_order, const MeshLib::Mesh& bulk_mesh); - -} // namespace ProcessLib diff --git a/ProcessLib/BoundaryCondition/NormalTractionBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/NormalTractionBoundaryConditionLocalAssembler.h index 7f1315dced38b34461fd55ae631de1b4d8e2687d..368a45ef2e88a1f4ac8025411eb383f054109a61 100644 --- a/ProcessLib/BoundaryCondition/NormalTractionBoundaryConditionLocalAssembler.h +++ b/ProcessLib/BoundaryCondition/NormalTractionBoundaryConditionLocalAssembler.h @@ -21,21 +21,19 @@ namespace ProcessLib { namespace NormalTractionBoundaryCondition { -template <typename ShapeMatricesTypeDisplacement, int GlobalDim, int NPoints> +template <typename ShapeMatricesType> struct IntegrationPointData final { IntegrationPointData( - typename ShapeMatricesTypeDisplacement::template VectorType< - NPoints * GlobalDim> const& Nu_times_n_, + typename ShapeMatricesType::ShapeMatrices::ShapeType const N_, + typename ShapeMatricesType::GlobalDimVectorType const n_, double const integration_weight_) - : Nu_times_n(Nu_times_n_), integration_weight(integration_weight_) + : N(N_), n(n_), integration_weight(integration_weight_) { } - // Shape matrix (for displacement) times element's normal vector. - typename ShapeMatricesTypeDisplacement::template VectorType< - NPoints * GlobalDim> const Nu_times_n; - + typename ShapeMatricesType::ShapeMatrices::ShapeType const N; + typename ShapeMatricesType::GlobalDimVectorType const n; double const integration_weight; }; @@ -56,11 +54,10 @@ class NormalTractionBoundaryConditionLocalAssembler final : public NormalTractionBoundaryConditionLocalAssemblerInterface { public: - using ShapeMatricesTypeDisplacement = + using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunctionDisplacement, GlobalDim>; - using GlobalDimVectorType = - typename ShapeMatrixPolicyType<ShapeFunctionDisplacement, - GlobalDim>::GlobalDimVectorType; + using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType; + using NodalVectorType = typename ShapeMatricesType::NodalVectorType; NormalTractionBoundaryConditionLocalAssembler( MeshLib::Element const& e, @@ -68,7 +65,9 @@ public: bool const is_axially_symmetric, unsigned const integration_order, Parameter<double> const& pressure) - : _integration_method(integration_order), _pressure(pressure) + : _integration_method(integration_order), + _pressure(pressure), + _element(e) { _local_rhs.setZero(local_matrix_size); @@ -78,10 +77,9 @@ public: _ip_data.reserve(n_integration_points); auto const shape_matrices_u = - initShapeMatrices<ShapeFunctionDisplacement, - ShapeMatricesTypeDisplacement, IntegrationMethod, - GlobalDim>(e, is_axially_symmetric, - _integration_method); + initShapeMatrices<ShapeFunctionDisplacement, ShapeMatricesType, + IntegrationMethod, GlobalDim>( + e, is_axially_symmetric, _integration_method); GlobalDimVectorType element_normal(GlobalDim); @@ -104,24 +102,13 @@ public: for (unsigned ip = 0; ip < n_integration_points; ip++) { - typename ShapeMatricesTypeDisplacement::template MatrixType< - GlobalDim, displacement_size> - N_u = ShapeMatricesTypeDisplacement::template MatrixType< - GlobalDim, displacement_size>::Zero(GlobalDim, - displacement_size); - for (int i = 0; i < static_cast<int>(GlobalDim); ++i) - { - N_u.template block<1, displacement_size / GlobalDim>( - i, i * displacement_size / GlobalDim) - .noalias() = shape_matrices_u[ip].N; - } double const integration_weight = _integration_method.getWeightedPoint(ip).getWeight() * shape_matrices_u[ip].integralMeasure * shape_matrices_u[ip].detJ; - _ip_data.emplace_back(N_u.transpose() * element_normal, + _ip_data.emplace_back(shape_matrices_u[ip].N, element_normal, integration_weight); } } @@ -137,18 +124,28 @@ public: unsigned const n_integration_points = _integration_method.getNumberOfPoints(); - SpatialPosition pos; - pos.setElementID(id); + NodalVectorType pressure = + _pressure.getNodalValuesOnElement(_element, t); for (unsigned ip = 0; ip < n_integration_points; ip++) { - pos.setIntegrationPoint(ip); - auto const& w = _ip_data[ip].integration_weight; - auto const& Nu_times_n = _ip_data[ip].Nu_times_n; + auto const& N = _ip_data[ip].N; + auto const& n = _ip_data[ip].n; - _local_rhs.noalias() -= - Nu_times_n.transpose() * _pressure(t, pos)[0] * w; + typename ShapeMatricesType::template MatrixType<GlobalDim, + displacement_size> + N_u = ShapeMatricesType::template MatrixType< + GlobalDim, displacement_size>::Zero(GlobalDim, + displacement_size); + for (int i = 0; i < static_cast<int>(GlobalDim); ++i) + { + N_u.template block<1, displacement_size / GlobalDim>( + i, i * displacement_size / GlobalDim) + .noalias() = N; + } + + _local_rhs.noalias() -= n.transpose() * N_u * pressure.dot(N) * w; } auto const indices = NumLib::getIndices(id, dof_table_boundary); @@ -161,17 +158,16 @@ private: static const int displacement_size = ShapeFunctionDisplacement::NPOINTS * GlobalDim; - std::vector<IntegrationPointData<ShapeMatricesTypeDisplacement, GlobalDim, - ShapeFunctionDisplacement::NPOINTS>, - Eigen::aligned_allocator<IntegrationPointData< - ShapeMatricesTypeDisplacement, GlobalDim, - ShapeFunctionDisplacement::NPOINTS>>> + std::vector< + IntegrationPointData<ShapeMatricesType>, + Eigen::aligned_allocator<IntegrationPointData<ShapeMatricesType>>> _ip_data; - typename ShapeMatricesTypeDisplacement::template VectorType< - displacement_size> + typename ShapeMatricesType::template VectorType<displacement_size> _local_rhs; + MeshLib::Element const& _element; + public: EIGEN_MAKE_ALIGNED_OPERATOR_NEW; }; diff --git a/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionLocalAssembler.h index ca87bdcd3a626274cdbf7ce592c02e6c64f30457..40437fa3317ff8863344972b7846068a8a1ccee6 100644 --- a/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionLocalAssembler.h +++ b/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionLocalAssembler.h @@ -41,9 +41,7 @@ public: bool is_axially_symmetric, unsigned const integration_order, PythonBoundaryConditionData const& data) - : Base(e, is_axially_symmetric, integration_order), - _data(data), - _element(e) + : Base(e, is_axially_symmetric, integration_order), _data(data) { } @@ -57,12 +55,12 @@ public: using FemType = NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>; FemType fe(*static_cast<const typename ShapeFunction::MeshElement*>( - &_element)); + &this->_element)); unsigned const num_integration_points = Base::_integration_method.getNumberOfPoints(); auto const num_var = _data.dof_table_bulk.getNumberOfVariables(); - auto const num_nodes = _element.getNumberOfNodes(); + auto const num_nodes = Base::_element.getNumberOfNodes(); auto const num_comp_total = _data.dof_table_bulk.getNumberOfComponents(); @@ -85,7 +83,8 @@ public: for (unsigned element_node_id = 0; element_node_id < num_nodes; ++element_node_id) { - auto const* const node = _element.getNode(element_node_id); + auto const* const node = + Base::_element.getNode(element_node_id); auto const boundary_node_id = node->getID(); auto const bulk_node_id = bulk_node_ids_map[boundary_node_id]; @@ -119,12 +118,12 @@ public: for (unsigned ip = 0; ip < num_integration_points; ip++) { - auto const& sm = Base::_shape_matrices[ip]; - auto const coords = fe.interpolateCoordinates(sm.N); + auto const& N = Base::_ns_and_weights[ip].N; + auto const& w = Base::_ns_and_weights[ip].weight; + auto const coords = fe.interpolateCoordinates(N); prim_vars = - sm.N * - primary_variables_mat; // Assumption: all primary variables - // have same shape functions. + N * primary_variables_mat; // Assumption: all primary variables + // have same shape functions. auto const flag_flux_dFlux = _data.bc_object->getFlux(t, coords, prim_vars_data); if (!_data.bc_object->isOverriddenNatural()) @@ -143,9 +142,7 @@ public: auto const flux = std::get<1>(flag_flux_dFlux); auto const& dFlux = std::get<2>(flag_flux_dFlux); - auto const& wp = Base::_integration_method.getWeightedPoint(ip); - auto const w = sm.detJ * wp.getWeight() * sm.integralMeasure; - local_rhs.noalias() += sm.N * (flux * w); + local_rhs.noalias() += N * (flux * w); if (static_cast<int>(dFlux.size()) != num_comp_total) { @@ -170,7 +167,7 @@ public: // The assignement -= takes into account the sign convention // of 1st-order in time ODE systems in OpenGeoSys. local_Jac.block(top, left, width, height).noalias() -= - sm.N.transpose() * (dFlux[comp] * w) * sm.N; + N.transpose() * (dFlux[comp] * w) * N; } } } @@ -193,7 +190,6 @@ public: private: PythonBoundaryConditionData const& _data; - MeshLib::Element const& _element; }; } // namespace ProcessLib diff --git a/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h index 715cbc137ca4a7d4cdedf35010c04e1d9e0a3214..434afb14b34402374a9f42c87d7d92cb350d0779 100644 --- a/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h +++ b/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h @@ -54,25 +54,22 @@ public: unsigned const n_integration_points = Base::_integration_method.getNumberOfPoints(); - SpatialPosition pos; - pos.setElementID(id); + typename Base::NodalVectorType const alpha = + _data.alpha.getNodalValuesOnElement(Base::_element, t); + typename Base::NodalVectorType const u_0 = + _data.u_0.getNodalValuesOnElement(Base::_element, t); for (unsigned ip = 0; ip < n_integration_points; ++ip) { - pos.setIntegrationPoint(ip); - auto const& sm = Base::_shape_matrices[ip]; - auto const& wp = Base::_integration_method.getWeightedPoint(ip); - - double const alpha = _data.alpha(t, pos)[0]; - double const u_0 = _data.u_0(t, pos)[0]; + auto const& ip_data = Base::_ns_and_weights[ip]; + auto const& N = ip_data.N; + auto const& w = ip_data.weight; // 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 * alpha * sm.detJ * wp.getWeight() * sm.integralMeasure; - _local_rhs.noalias() += sm.N * alpha * u_0 * sm.detJ * - wp.getWeight() * sm.integralMeasure; + _local_K.diagonal().noalias() += N * alpha.dot(N) * w; + _local_rhs.noalias() += N * alpha.dot(N) * u_0.dot(N) * w; } auto const indices = NumLib::getIndices(id, dof_table_boundary); diff --git a/ProcessLib/BoundaryCondition/NonuniformVariableDependentNeumannBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/VariableDependentNeumannBoundaryCondition.cpp similarity index 51% rename from ProcessLib/BoundaryCondition/NonuniformVariableDependentNeumannBoundaryCondition.cpp rename to ProcessLib/BoundaryCondition/VariableDependentNeumannBoundaryCondition.cpp index 2af6cb375039f0f1f4710a61976e631798936db3..36f30cd6753d1f4b4eb6fdeb092b11cd4dacbf9c 100644 --- a/ProcessLib/BoundaryCondition/NonuniformVariableDependentNeumannBoundaryCondition.cpp +++ b/ProcessLib/BoundaryCondition/VariableDependentNeumannBoundaryCondition.cpp @@ -7,69 +7,59 @@ * */ -#include "NonuniformVariableDependentNeumannBoundaryCondition.h" +#include "VariableDependentNeumannBoundaryCondition.h" -#include "MeshLib/IO/readMeshFromFile.h" #include "ProcessLib/Utils/ProcessUtils.h" namespace ProcessLib { -std::unique_ptr<NonuniformVariableDependentNeumannBoundaryCondition> -createNonuniformVariableDependentNeumannBoundaryCondition( - BaseLib::ConfigTree const& config, MeshLib::Mesh const& boundary_mesh, +std::unique_ptr<VariableDependentNeumannBoundaryCondition> +createVariableDependentNeumannBoundaryCondition( + BaseLib::ConfigTree const& config, MeshLib::Mesh const& bc_mesh, NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id, int const component_id, unsigned const integration_order, - unsigned const shapefunction_order, MeshLib::Mesh const& bulk_mesh) + unsigned const shapefunction_order, unsigned const global_dim, + std::vector<std::unique_ptr<ParameterBase>> const& parameters) { - DBUG("Constructing NonuniformVariableDependentNeumann BC from config."); + DBUG("Constructing VariableDependentNeumann BC from config."); //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__type} - config.checkConfigParameter("type", "NonuniformVariableDependentNeumann"); + config.checkConfigParameter("type", "VariableDependentNeumann"); if (dof_table.getNumberOfVariables() != 2) { OGS_FATAL( - "NonuniformVariableDependentNeumann BC only implemented for 2 " + "VariableDependentNeumann BC only implemented for 2 " "variable processes."); } assert(variable_id == 0 || variable_id == 1); - // TODO finally use ProcessLib::Parameter here auto const constant_name = - //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__NonuniformVariableDependentNeumann__constant_name} + //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__VariableDependentNeumann__constant_name} config.getConfigParameter<std::string>("constant_name"); - auto const& constant = - *boundary_mesh.getProperties().getPropertyVector<double>( - constant_name, MeshLib::MeshItemType::Node, 1); + auto const& constant = findParameter<double>(constant_name, parameters, 1); auto const coefficient_current_variable_name = - //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__NonuniformVariableDependentNeumann__coefficient_current_variable_name} + //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__VariableDependentNeumann__coefficient_current_variable_name} config.getConfigParameter<std::string>( "coefficient_current_variable_name"); auto const& coefficient_current_variable = - *boundary_mesh.getProperties().getPropertyVector<double>( - coefficient_current_variable_name, MeshLib::MeshItemType::Node, 1); + findParameter<double>(coefficient_current_variable_name, parameters, 1); auto const coefficient_other_variable_name = - //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__NonuniformVariableDependentNeumann__coefficient_other_variable_name} + //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__VariableDependentNeumann__coefficient_other_variable_name} config.getConfigParameter<std::string>( "coefficient_other_variable_name"); auto const& coefficient_other_variable = - *boundary_mesh.getProperties().getPropertyVector<double>( - coefficient_other_variable_name, MeshLib::MeshItemType::Node, 1); + findParameter<double>(coefficient_other_variable_name, parameters, 1); auto const coefficient_mixed_variables_name = - //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__NonuniformVariableDependentNeumann__coefficient_mixed_variables_name} + //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__VariableDependentNeumann__coefficient_mixed_variables_name} config.getConfigParameter<std::string>( "coefficient_mixed_variables_name"); auto const& coefficient_mixed_variables = - *boundary_mesh.getProperties().getPropertyVector<double>( - coefficient_mixed_variables_name, MeshLib::MeshItemType::Node, 1); + findParameter<double>(coefficient_mixed_variables_name, parameters, 1); - std::string const mapping_to_bulk_nodes_property = "bulk_node_ids"; - boundary_mesh.getProperties().getPropertyVector<std::size_t>( - mapping_to_bulk_nodes_property, MeshLib::MeshItemType::Node, 1); - - std::vector<MeshLib::Node*> const& bc_nodes = boundary_mesh.getNodes(); - MeshLib::MeshSubset bc_mesh_subset(boundary_mesh, bc_nodes); + std::vector<MeshLib::Node*> const& bc_nodes = bc_mesh.getNodes(); + MeshLib::MeshSubset bc_mesh_subset(bc_mesh, bc_nodes); auto const& dof_table_boundary_other_variable = *dof_table.deriveBoundaryConstrainedMap( (variable_id + 1) % 2, {component_id}, std::move(bc_mesh_subset)); @@ -81,19 +71,17 @@ createNonuniformVariableDependentNeumannBoundaryCondition( // parameters are not read and will cause an error. // TODO (naumov): Add a function to ConfigTree for skipping the tags of the // subtree and move the code up in createBoundaryCondition(). - if (boundary_mesh.getDimension() == 0 && - boundary_mesh.getNumberOfNodes() == 0 && - boundary_mesh.getNumberOfElements() == 0) + if (bc_mesh.getDimension() == 0 && bc_mesh.getNumberOfNodes() == 0 && + bc_mesh.getNumberOfElements() == 0) { return nullptr; } #endif // USE_PETSC - return std::make_unique< - NonuniformVariableDependentNeumannBoundaryCondition>( + return std::make_unique<VariableDependentNeumannBoundaryCondition>( integration_order, shapefunction_order, dof_table, variable_id, - component_id, bulk_mesh.getDimension(), boundary_mesh, - NonuniformVariableDependentNeumannBoundaryConditionData{ + component_id, global_dim, bc_mesh, + VariableDependentNeumannBoundaryConditionData{ constant, coefficient_current_variable, coefficient_other_variable, coefficient_mixed_variables, dof_table_boundary_other_variable}); } diff --git a/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.h b/ProcessLib/BoundaryCondition/VariableDependentNeumannBoundaryCondition.h similarity index 50% rename from ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.h rename to ProcessLib/BoundaryCondition/VariableDependentNeumannBoundaryCondition.h index 50cb24c44bbb2fecef3e63215f83fc149e864a40..82f1170ff5c824b5cad03e4ddfcabb119c3e8e42 100644 --- a/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.h +++ b/ProcessLib/BoundaryCondition/VariableDependentNeumannBoundaryCondition.h @@ -9,22 +9,23 @@ #pragma once -#include "GenericNonuniformNaturalBoundaryCondition.h" +#include "GenericNaturalBoundaryCondition.h" #include "MeshLib/PropertyVector.h" -#include "NonuniformNeumannBoundaryConditionLocalAssembler.h" +#include "VariableDependentNeumannBoundaryConditionLocalAssembler.h" namespace ProcessLib { -using NonuniformNeumannBoundaryCondition = - GenericNonuniformNaturalBoundaryCondition< - NonuniformNeumannBoundaryConditionData, - NonuniformNeumannBoundaryConditionLocalAssembler>; +using VariableDependentNeumannBoundaryCondition = + GenericNaturalBoundaryCondition< + VariableDependentNeumannBoundaryConditionData, + VariableDependentNeumannBoundaryConditionLocalAssembler>; -std::unique_ptr<NonuniformNeumannBoundaryCondition> -createNonuniformNeumannBoundaryCondition( +std::unique_ptr<VariableDependentNeumannBoundaryCondition> +createVariableDependentNeumannBoundaryCondition( BaseLib::ConfigTree const& config, MeshLib::Mesh const& boundary_mesh, NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id, int const component_id, unsigned const integration_order, - unsigned const shapefunction_order, const MeshLib::Mesh& bulk_mesh); + unsigned const shapefunction_order, unsigned const global_dim, + std::vector<std::unique_ptr<ParameterBase>> const& parameters); } // namespace ProcessLib diff --git a/ProcessLib/BoundaryCondition/NonuniformVariableDependentNeumannBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/VariableDependentNeumannBoundaryConditionLocalAssembler.h similarity index 68% rename from ProcessLib/BoundaryCondition/NonuniformVariableDependentNeumannBoundaryConditionLocalAssembler.h rename to ProcessLib/BoundaryCondition/VariableDependentNeumannBoundaryConditionLocalAssembler.h index df2972c85c5bcde7225c48a2c01bf62428363dbb..b5e0b0112f1c4c736d1796e82f0e301e1e46cc16 100644 --- a/ProcessLib/BoundaryCondition/NonuniformVariableDependentNeumannBoundaryConditionLocalAssembler.h +++ b/ProcessLib/BoundaryCondition/VariableDependentNeumannBoundaryConditionLocalAssembler.h @@ -14,28 +14,27 @@ #include "NumLib/Function/Interpolation.h" #include "ProcessLib/Parameter/MeshNodeParameter.h" -#include "GenericNonuniformNaturalBoundaryConditionLocalAssembler.h" +#include "GenericNaturalBoundaryConditionLocalAssembler.h" namespace ProcessLib { -struct NonuniformVariableDependentNeumannBoundaryConditionData +struct VariableDependentNeumannBoundaryConditionData { - /* TODO use Parameter */ - MeshLib::PropertyVector<double> const& constant; - MeshLib::PropertyVector<double> const& coefficient_current_variable; - MeshLib::PropertyVector<double> const& coefficient_other_variable; - MeshLib::PropertyVector<double> const& coefficient_mixed_variables; + Parameter<double> const& constant; + Parameter<double> const& coefficient_current_variable; + Parameter<double> const& coefficient_other_variable; + Parameter<double> const& coefficient_mixed_variables; // Used for mapping boundary nodes to bulk nodes. NumLib::LocalToGlobalIndexMap const& dof_table_boundary_other_variable; }; template <typename ShapeFunction, typename IntegrationMethod, unsigned GlobalDim> -class NonuniformVariableDependentNeumannBoundaryConditionLocalAssembler final - : public GenericNonuniformNaturalBoundaryConditionLocalAssembler< +class VariableDependentNeumannBoundaryConditionLocalAssembler final + : public GenericNaturalBoundaryConditionLocalAssembler< ShapeFunction, IntegrationMethod, GlobalDim> { - using Base = GenericNonuniformNaturalBoundaryConditionLocalAssembler< + using Base = GenericNaturalBoundaryConditionLocalAssembler< ShapeFunction, IntegrationMethod, GlobalDim>; using NodalVectorType = typename Base::NodalVectorType; using NodalMatrixType = typename Base::NodalMatrixType; @@ -43,12 +42,12 @@ class NonuniformVariableDependentNeumannBoundaryConditionLocalAssembler final public: /// The neumann_bc_term factor is directly integrated into the local /// element matrix. - NonuniformVariableDependentNeumannBoundaryConditionLocalAssembler( + VariableDependentNeumannBoundaryConditionLocalAssembler( MeshLib::Element const& e, std::size_t const local_matrix_size, bool const is_axially_symmetric, unsigned const integration_order, - NonuniformVariableDependentNeumannBoundaryConditionData const& data) + VariableDependentNeumannBoundaryConditionData const& data) : Base(e, is_axially_symmetric, integration_order), _data(data), _local_matrix_size(local_matrix_size) @@ -62,28 +61,21 @@ public: { NodalVectorType _local_rhs(_local_matrix_size); _local_rhs.setZero(); - MeshNodeParameter<double> constant_values{"ConstantValues", - _data.constant}; - MeshNodeParameter<double> coefficient_current_variable_values{ - "CurrentVariableValues", _data.coefficient_current_variable}; - MeshNodeParameter<double> coefficient_other_variable_values{ - "OtherVariableValues", _data.coefficient_other_variable}; - MeshNodeParameter<double> coefficient_mixed_variables_values{ - "MixedVariablesValues", _data.coefficient_mixed_variables}; // Get element nodes for the interpolation from nodes to // integration point. NodalVectorType const constant_node_values = - constant_values.getNodalValuesOnElement(Base::_element, t); + _data.constant.getNodalValuesOnElement(Base::_element, t); NodalVectorType const coefficient_current_variable_node_values = - coefficient_current_variable_values.getNodalValuesOnElement( + _data.coefficient_current_variable.getNodalValuesOnElement( Base::_element, t); NodalVectorType const coefficient_other_variable_node_values = - coefficient_other_variable_values.getNodalValuesOnElement( + _data.coefficient_other_variable.getNodalValuesOnElement( Base::_element, t); NodalVectorType const coefficient_mixed_variables_node_values = - coefficient_mixed_variables_values.getNodalValuesOnElement( + _data.coefficient_mixed_variables.getNodalValuesOnElement( Base::_element, t); - unsigned const n_integration_points = _local_matrix_size; + unsigned const n_integration_points = + Base::_integration_method.getNumberOfPoints(); auto const indices_current_variable = NumLib::getIndices(mesh_item_id, dof_table_boundary); @@ -121,7 +113,7 @@ public: } private: - NonuniformVariableDependentNeumannBoundaryConditionData const& _data; + VariableDependentNeumannBoundaryConditionData const& _data; std::size_t const _local_matrix_size; }; diff --git a/ProcessLib/Parameter/CurveScaledParameter.h b/ProcessLib/Parameter/CurveScaledParameter.h index 3e4babb03045cd949e84faf603e6b08d1fe71a72..bce7cfbedd3c7dfe38c8bf955d25e74e93e1d3ac 100644 --- a/ProcessLib/Parameter/CurveScaledParameter.h +++ b/ProcessLib/Parameter/CurveScaledParameter.h @@ -35,6 +35,7 @@ struct CurveScaledParameter final : public Parameter<T> { { _parameter = &findParameter<T>(_referenced_parameter_name, parameters, 0); + ParameterBase::_mesh = _parameter->mesh(); } int getNumberOfComponents() const override diff --git a/ProcessLib/Parameter/FunctionParameter.h b/ProcessLib/Parameter/FunctionParameter.h index 82b84338d58f7148e239bf91dc610e3d99968928..f32a3cddc236343f82ad888904d9365f217c414a 100644 --- a/ProcessLib/Parameter/FunctionParameter.h +++ b/ProcessLib/Parameter/FunctionParameter.h @@ -39,15 +39,15 @@ struct FunctionParameter final : public Parameter<T> /** * Constructing from a vector of expressions * - * @param name_ the parameter's name - * @param mesh_ a mesh object - * @param vec_expression_str_ a vector of mathematical expressions + * @param name the parameter's name + * @param mesh the parameter's domain of definition. + * @param vec_expression_str a vector of mathematical expressions * The vector size specifies the number of components of the parameter. */ - FunctionParameter(std::string const& name_, - MeshLib::Mesh const& mesh_, - std::vector<std::string> const& vec_expression_str_) - : Parameter<T>(name_), _mesh(mesh_), _vec_expression_str(vec_expression_str_) + FunctionParameter(std::string const& name, + MeshLib::Mesh const& mesh, + std::vector<std::string> const& vec_expression_str) + : Parameter<T>(name, &mesh), _vec_expression_str(vec_expression_str) { _symbol_table.add_constants(); _symbol_table.create_variable("x"); @@ -90,7 +90,8 @@ struct FunctionParameter final : public Parameter<T> } else if (pos.getNodeID()) { - auto const& node = *_mesh.getNode(pos.getNodeID().get()); + auto const& node = + *ParameterBase::_mesh->getNode(pos.getNodeID().get()); x = node[0]; y = node[1]; z = node[2]; @@ -105,7 +106,6 @@ struct FunctionParameter final : public Parameter<T> } private: - MeshLib::Mesh const& _mesh; std::vector<std::string> const _vec_expression_str; symbol_table_t _symbol_table; std::vector<expression_t> _vec_expression; diff --git a/ProcessLib/Parameter/GroupBasedParameter.cpp b/ProcessLib/Parameter/GroupBasedParameter.cpp index f6b256b233530dc06c5439f4af8291fb023d10fe..e939846c857bb92794305146a3be0f224026f23f 100644 --- a/ProcessLib/Parameter/GroupBasedParameter.cpp +++ b/ProcessLib/Parameter/GroupBasedParameter.cpp @@ -99,13 +99,13 @@ std::unique_ptr<ParameterBase> createGroupBasedParameter( { return std::make_unique< GroupBasedParameter<double, MeshLib::MeshItemType::Node>>( - name, *group_id_property, vec_values); + name, mesh, *group_id_property, vec_values); } if (group_id_property->getMeshItemType() == MeshLib::MeshItemType::Cell) { return std::make_unique< GroupBasedParameter<double, MeshLib::MeshItemType::Cell>>( - name, *group_id_property, vec_values); + name, mesh, *group_id_property, vec_values); } OGS_FATAL("Mesh item type of the specified property is not supported."); diff --git a/ProcessLib/Parameter/GroupBasedParameter.h b/ProcessLib/Parameter/GroupBasedParameter.h index 67b7ecfcf6fead611e5dfd18bb42f3cfeaba7ece..1a1d217c909c3ce5bf96246eb944179ab87c6540 100644 --- a/ProcessLib/Parameter/GroupBasedParameter.h +++ b/ProcessLib/Parameter/GroupBasedParameter.h @@ -36,14 +36,16 @@ struct GroupBasedParameter final * Constructing from a property vector of index and corresponding values * * @param name_ the parameter's name + * @param mesh the parameter's domain of definition. * @param property a property vector of index for mesh items * @param vec_values a vector of values for each index */ GroupBasedParameter(std::string const& name_, + MeshLib::Mesh const& mesh, MeshLib::PropertyVector<int> const& property, std::vector<std::vector<double>> vec_values) - : Parameter<T>(name_), + : Parameter<T>(name_, &mesh), _property_index(property), _vec_values(std::move(vec_values)) { diff --git a/ProcessLib/Parameter/MeshElementParameter.cpp b/ProcessLib/Parameter/MeshElementParameter.cpp index aaa20760daabe2547be608b521c9b409c9edf835..8c707397d762ec99737dbcc916871b4eabe2e35f 100644 --- a/ProcessLib/Parameter/MeshElementParameter.cpp +++ b/ProcessLib/Parameter/MeshElementParameter.cpp @@ -32,7 +32,8 @@ std::unique_ptr<ParameterBase> createMeshElementParameter( field_name.c_str()); } - return std::make_unique<MeshElementParameter<double>>(name, *property); + return std::make_unique<MeshElementParameter<double>>(name, mesh, + *property); } } // namespace ProcessLib diff --git a/ProcessLib/Parameter/MeshElementParameter.h b/ProcessLib/Parameter/MeshElementParameter.h index 1146105622d4fe12c8977ae5970bb94ef338d097..1999f6f014b7eb563df6a2f80bebe72c69975f91 100644 --- a/ProcessLib/Parameter/MeshElementParameter.h +++ b/ProcessLib/Parameter/MeshElementParameter.h @@ -23,9 +23,9 @@ namespace ProcessLib template <typename T> struct MeshElementParameter final : public Parameter<T> { MeshElementParameter(std::string const& name_, + MeshLib::Mesh const& mesh, MeshLib::PropertyVector<T> const& property) - : Parameter<T>(name_), - _property(property) + : Parameter<T>(name_, &mesh), _property(property) { } diff --git a/ProcessLib/Parameter/MeshNodeParameter.cpp b/ProcessLib/Parameter/MeshNodeParameter.cpp index 48fc8f5a281de41edd95e91a6145056f32da9d3e..11b9ecbabe2186d6c82605503bf7ad8560c262ea 100644 --- a/ProcessLib/Parameter/MeshNodeParameter.cpp +++ b/ProcessLib/Parameter/MeshNodeParameter.cpp @@ -32,7 +32,7 @@ std::unique_ptr<ParameterBase> createMeshNodeParameter( field_name.c_str()); } - return std::make_unique<MeshNodeParameter<double>>(name, *property); + return std::make_unique<MeshNodeParameter<double>>(name, mesh, *property); } } // namespace ProcessLib diff --git a/ProcessLib/Parameter/MeshNodeParameter.h b/ProcessLib/Parameter/MeshNodeParameter.h index 5c013c22dcf87e60395b2489022091d78c3d2aea..5ea4da35a38b86bfbb3924d83d1c3459d2521a0f 100644 --- a/ProcessLib/Parameter/MeshNodeParameter.h +++ b/ProcessLib/Parameter/MeshNodeParameter.h @@ -27,9 +27,9 @@ namespace ProcessLib template <typename T> struct MeshNodeParameter final : public Parameter<T> { MeshNodeParameter(std::string const& name_, + MeshLib::Mesh const& mesh, MeshLib::PropertyVector<T> const& property) - : Parameter<T>(name_), - _property(property) + : Parameter<T>(name_, &mesh), _property(property) { } diff --git a/ProcessLib/Parameter/Parameter.cpp b/ProcessLib/Parameter/Parameter.cpp index 8d21a617f5c826b559ec5932ecacb0f8e8648c03..2cddcff1c57459e195d3c1b417d94cf0e8d25c9d 100644 --- a/ProcessLib/Parameter/Parameter.cpp +++ b/ProcessLib/Parameter/Parameter.cpp @@ -33,6 +33,17 @@ std::unique_ptr<ParameterBase> createParameter( //! \ogs_file_param{prj__parameters__parameter__type} auto const type = config.peekConfigParameter<std::string>("type"); + // Either the mesh name is given, or the first mesh's name will be + // taken. + auto const mesh_name = + //! \ogs_file_param{prj__parameters__parameter__mesh} + config.getConfigParameter<std::string>("mesh", meshes[0]->getName()); + + auto const& mesh = *BaseLib::findElementOrError( + begin(meshes), end(meshes), + [&mesh_name](auto const& m) { return m->getName() == mesh_name; }, + "Expected to find a mesh named " + mesh_name + "."); + // Create parameter based on the provided type. if (type == "Constant") { @@ -49,25 +60,25 @@ std::unique_ptr<ParameterBase> createParameter( if (type == "Function") { INFO("FunctionParameter: %s", name.c_str()); - auto param = createFunctionParameter(name, config, *meshes.front()); + auto param = createFunctionParameter(name, config, mesh); return param; } if (type == "Group") { INFO("GroupBasedParameter: %s", name.c_str()); - auto param = createGroupBasedParameter(name, config, *meshes.front()); + auto param = createGroupBasedParameter(name, config, mesh); return param; } if (type == "MeshElement") { INFO("MeshElementParameter: %s", name.c_str()); - auto param = createMeshElementParameter(name, config, *meshes.front()); + auto param = createMeshElementParameter(name, config, mesh); return param; } if (type == "MeshNode") { INFO("MeshNodeParameter: %s", name.c_str()); - auto param = createMeshNodeParameter(name, config, *meshes.front()); + auto param = createMeshNodeParameter(name, config, mesh); return param; } diff --git a/ProcessLib/Parameter/Parameter.h b/ProcessLib/Parameter/Parameter.h index e350d341550b7d21ea9d4e68572b2e864549374f..51c71a82d2f3e888713827ca952b80aa39ea912a 100644 --- a/ProcessLib/Parameter/Parameter.h +++ b/ProcessLib/Parameter/Parameter.h @@ -44,7 +44,12 @@ namespace ProcessLib /// Its property name helps addressing the right parameter. struct ParameterBase { - ParameterBase(std::string name_) : name(std::move(name_)) {} + explicit ParameterBase(std::string name_, + MeshLib::Mesh const* mesh = nullptr) + : name(std::move(name_)), _mesh(mesh) + { + } + virtual ~ParameterBase() = default; virtual bool isTimeDependent() const = 0; @@ -57,7 +62,14 @@ struct ParameterBase { } + MeshLib::Mesh const* mesh() const { return _mesh; } + std::string const name; + +protected: + /// A mesh on which the parameter is defined. Some parameters might be + /// mesh-independent. + MeshLib::Mesh const* _mesh; }; /*! A Parameter is a function \f$ (t, x) \mapsto f(t, x) \in T^n \f$. @@ -69,7 +81,8 @@ struct ParameterBase template <typename T> struct Parameter : public ParameterBase { - Parameter(std::string const& name_) : ParameterBase(name_) {} + using ParameterBase::ParameterBase; + ~Parameter() override = default; //! Returns the number of components this Parameter has at every position and diff --git a/Tests/Data/Elliptic/nonuniform_bc_Groundwaterflow/dirichlet_nonuniform.prj b/Tests/Data/Elliptic/nonuniform_bc_Groundwaterflow/dirichlet_nonuniform.prj index 99c7e8e777e3e966e8267a66f431ed131a35336d..7518b9483bbb621aee9d51739beebd26cb73eab1 100644 --- a/Tests/Data/Elliptic/nonuniform_bc_Groundwaterflow/dirichlet_nonuniform.prj +++ b/Tests/Data/Elliptic/nonuniform_bc_Groundwaterflow/dirichlet_nonuniform.prj @@ -59,6 +59,12 @@ <type>Constant</type> <value>-1</value> </parameter> + <parameter> + <name>p_Dirichlet_left</name> + <type>MeshNode</type> + <mesh>square_1x1_left</mesh> + <field_name>linearY</field_name> + </parameter> </parameters> <process_variables> <process_variable> @@ -69,9 +75,9 @@ <initial_condition>p0</initial_condition> <boundary_conditions> <boundary_condition> - <type>NonuniformDirichlet</type> <mesh>square_1x1_left</mesh> - <field_name>linearY</field_name> + <type>Dirichlet</type> + <parameter>p_Dirichlet_left</parameter> </boundary_condition> <boundary_condition> <mesh>square_1x1_right</mesh> diff --git a/Tests/Data/Elliptic/nonuniform_bc_Groundwaterflow/inhomogeneous_permeability.prj b/Tests/Data/Elliptic/nonuniform_bc_Groundwaterflow/inhomogeneous_permeability.prj index 678a1e2bc9fc02b24d57ff70ec873d2d57fb3675..c30ee66bbe1267382f9e1e07fb3463b43e141d6d 100644 --- a/Tests/Data/Elliptic/nonuniform_bc_Groundwaterflow/inhomogeneous_permeability.prj +++ b/Tests/Data/Elliptic/nonuniform_bc_Groundwaterflow/inhomogeneous_permeability.prj @@ -61,6 +61,12 @@ <type>Constant</type> <value>0</value> </parameter> + <parameter> + <name>mass_flux</name> + <type>MeshNode</type> + <mesh>inhomogeneous_permeability_top</mesh> + <field_name>mass_flux</field_name> + </parameter> </parameters> <process_variables> <process_variable> @@ -72,8 +78,8 @@ <boundary_conditions> <boundary_condition> <mesh>inhomogeneous_permeability_top</mesh> - <type>NonuniformNeumann</type> - <field_name>mass_flux</field_name> + <type>Neumann</type> + <parameter>mass_flux</parameter> </boundary_condition> <boundary_condition> <mesh>inhomogeneous_permeability_bottom</mesh> diff --git a/Tests/Data/Elliptic/nonuniform_bc_Groundwaterflow/neumann_nonuniform.prj b/Tests/Data/Elliptic/nonuniform_bc_Groundwaterflow/neumann_nonuniform.prj index b64e5f1bcebebe056e2e2d202ad5468021caa404..b2e30d62990b93ff8f4649a89338a506eec143a9 100644 --- a/Tests/Data/Elliptic/nonuniform_bc_Groundwaterflow/neumann_nonuniform.prj +++ b/Tests/Data/Elliptic/nonuniform_bc_Groundwaterflow/neumann_nonuniform.prj @@ -61,6 +61,12 @@ <type>Constant</type> <value>-1</value> </parameter> + <parameter> + <name>p_Neumann_left</name> + <type>MeshNode</type> + <mesh>square_1x1_left</mesh> + <field_name>cosY</field_name> + </parameter> </parameters> <process_variables> <process_variable> @@ -71,9 +77,9 @@ <initial_condition>p0</initial_condition> <boundary_conditions> <boundary_condition> - <type>NonuniformNeumann</type> + <type>Neumann</type> <mesh>square_1x1_left</mesh> - <field_name>cosY</field_name> + <parameter>p_Neumann_left</parameter> </boundary_condition> <boundary_condition> <type>Dirichlet</type> diff --git a/Tests/Data/Parabolic/ComponentTransport/VariableNeumannBoundary/vdbc_input.prj b/Tests/Data/Parabolic/ComponentTransport/VariableNeumannBoundary/vdbc_input.prj index 46fed65132efb7234954f831e6df2b798610f1b6..8751123e5979b0e9db4f81bcb6a9b01048a02f77 100644 --- a/Tests/Data/Parabolic/ComponentTransport/VariableNeumannBoundary/vdbc_input.prj +++ b/Tests/Data/Parabolic/ComponentTransport/VariableNeumannBoundary/vdbc_input.prj @@ -169,6 +169,18 @@ <type>MeshNode</type> <field_name>p_ini</field_name> </parameter> + <parameter> + <name>c_left</name> + <mesh>vdbc_input_leftBoundary</mesh> + <type>MeshNode</type> + <field_name>c_ini</field_name> + </parameter> + <parameter> + <name>p_left</name> + <mesh>vdbc_input_leftBoundary</mesh> + <type>MeshNode</type> + <field_name>p_ini</field_name> + </parameter> <parameter> <name>constant_porosity_parameter</name> <type>Constant</type> @@ -179,6 +191,54 @@ <type>Constant</type> <values>1.e-9 0 0 1.e-9</values> </parameter> + <parameter> + <name>Cconstant</name> + <mesh>vdbc_input_rightBoundary</mesh> + <type>MeshNode</type> + <field_name>Cconstant</field_name> + </parameter> + <parameter> + <name>Cprefac1</name> + <mesh>vdbc_input_rightBoundary</mesh> + <type>MeshNode</type> + <field_name>Cprefac1</field_name> + </parameter> + <parameter> + <name>Cprefac2</name> + <mesh>vdbc_input_rightBoundary</mesh> + <type>MeshNode</type> + <field_name>Cprefac2</field_name> + </parameter> + <parameter> + <name>Cprefac3</name> + <mesh>vdbc_input_rightBoundary</mesh> + <type>MeshNode</type> + <field_name>Cprefac3</field_name> + </parameter> + <parameter> + <name>Pconstant</name> + <mesh>vdbc_input_rightBoundary</mesh> + <type>MeshNode</type> + <field_name>Pconstant</field_name> + </parameter> + <parameter> + <name>Pprefac1</name> + <mesh>vdbc_input_rightBoundary</mesh> + <type>MeshNode</type> + <field_name>Pprefac1</field_name> + </parameter> + <parameter> + <name>Pprefac2</name> + <mesh>vdbc_input_rightBoundary</mesh> + <type>MeshNode</type> + <field_name>Pprefac2</field_name> + </parameter> + <parameter> + <name>Pprefac3</name> + <mesh>vdbc_input_rightBoundary</mesh> + <type>MeshNode</type> + <field_name>Pprefac3</field_name> + </parameter> </parameters> <process_variables> <process_variable> @@ -188,7 +248,7 @@ <initial_condition>c_ini</initial_condition> <boundary_conditions> <boundary_condition> - <type>NonuniformVariableDependentNeumann</type> + <type>VariableDependentNeumann</type> <mesh>vdbc_input_rightBoundary</mesh> <constant_name>Cconstant</constant_name> <coefficient_current_variable_name>Cprefac1</coefficient_current_variable_name> @@ -196,9 +256,9 @@ <coefficient_mixed_variables_name>Cprefac3</coefficient_mixed_variables_name> </boundary_condition> <boundary_condition> - <type>NonuniformDirichlet</type> + <type>Dirichlet</type> <mesh>vdbc_input_leftBoundary</mesh> - <field_name>c_ini</field_name> + <parameter>c_left</parameter> </boundary_condition> </boundary_conditions> </process_variable> @@ -209,7 +269,7 @@ <initial_condition>p_ini</initial_condition> <boundary_conditions> <boundary_condition> - <type>NonuniformVariableDependentNeumann</type> + <type>VariableDependentNeumann</type> <mesh>vdbc_input_rightBoundary</mesh> <constant_name>Pconstant</constant_name> <coefficient_current_variable_name>Pprefac1</coefficient_current_variable_name> @@ -217,9 +277,9 @@ <coefficient_mixed_variables_name>Pprefac3</coefficient_mixed_variables_name> </boundary_condition> <boundary_condition> - <type>NonuniformDirichlet</type> + <type>Dirichlet</type> <mesh>vdbc_input_leftBoundary</mesh> - <field_name>p_ini</field_name> + <parameter>p_left</parameter> </boundary_condition> </boundary_conditions> </process_variable> diff --git a/Tests/Data/Parabolic/ComponentTransport/goswami/goswami_input.prj b/Tests/Data/Parabolic/ComponentTransport/goswami/goswami_input.prj index 0609c824eeab40229b7963e2d2d044d5b86b75a6..54b45b696c879c06716b8be8212706d8a64384dc 100644 --- a/Tests/Data/Parabolic/ComponentTransport/goswami/goswami_input.prj +++ b/Tests/Data/Parabolic/ComponentTransport/goswami/goswami_input.prj @@ -181,6 +181,30 @@ <type>Constant</type> <values>1.2388e-9 0 0 1.2388e-9</values> </parameter> + <parameter> + <name>concentration_left</name> + <type>MeshNode</type> + <mesh>goswami_input_leftBoundary</mesh> + <field_name>c_ini</field_name> + </parameter> + <parameter> + <name>concentration_right</name> + <type>MeshNode</type> + <mesh>goswami_input_rightBoundary</mesh> + <field_name>c_ini</field_name> + </parameter> + <parameter> + <name>pressure_left</name> + <type>MeshNode</type> + <mesh>goswami_input_leftBoundary</mesh> + <field_name>p_ini</field_name> + </parameter> + <parameter> + <name>pressure_right</name> + <type>MeshNode</type> + <mesh>goswami_input_rightBoundary</mesh> + <field_name>p_ini</field_name> + </parameter> </parameters> <process_variables> <process_variable> @@ -190,14 +214,14 @@ <initial_condition>c_ini</initial_condition> <boundary_conditions> <boundary_condition> - <type>NonuniformDirichlet</type> + <type>Dirichlet</type> <mesh>goswami_input_leftBoundary</mesh> - <field_name>c_ini</field_name> + <parameter>concentration_left</parameter> </boundary_condition> <boundary_condition> - <type>NonuniformDirichlet</type> + <type>Dirichlet</type> <mesh>goswami_input_rightBoundary</mesh> - <field_name>c_ini</field_name> + <parameter>concentration_right</parameter> </boundary_condition> </boundary_conditions> </process_variable> @@ -208,14 +232,14 @@ <initial_condition>p_ini</initial_condition> <boundary_conditions> <boundary_condition> - <type>NonuniformDirichlet</type> + <type>Dirichlet</type> <mesh>goswami_input_leftBoundary</mesh> - <field_name>p_ini</field_name> + <parameter>pressure_left</parameter> </boundary_condition> <boundary_condition> - <type>NonuniformDirichlet</type> + <type>Dirichlet</type> <mesh>goswami_input_rightBoundary</mesh> - <field_name>p_ini</field_name> + <parameter>pressure_right</parameter> </boundary_condition> </boundary_conditions> </process_variable> diff --git a/Tests/Data/Parabolic/ComponentTransport/heterogeneous/ogs5_H_3D/ogs5_H_3d.prj b/Tests/Data/Parabolic/ComponentTransport/heterogeneous/ogs5_H_3D/ogs5_H_3d.prj index 2e84856f6cefaca1a00d8de681d94d44750783fd..6471307bac285d34db23ab67ff04988f7ce34785 100644 --- a/Tests/Data/Parabolic/ComponentTransport/heterogeneous/ogs5_H_3D/ogs5_H_3d.prj +++ b/Tests/Data/Parabolic/ComponentTransport/heterogeneous/ogs5_H_3D/ogs5_H_3d.prj @@ -181,6 +181,18 @@ <!--type>Constant</type> <values>1.239e-7 0 0 1.239e-7</values--> </parameter> + <parameter> + <name>bc_left</name> + <type>MeshNode</type> + <mesh>3D_HGW_1_l</mesh> + <field_name>p_bc</field_name> + </parameter> + <parameter> + <name>bc_right</name> + <type>MeshNode</type> + <mesh>3D_HGW_1_r</mesh> + <field_name>p_bc</field_name> + </parameter> </parameters> <process_variables> <process_variable> @@ -208,14 +220,14 @@ <initial_condition>p_ini</initial_condition> <boundary_conditions> <boundary_condition> - <type>NonuniformDirichlet</type> + <type>Dirichlet</type> <mesh>3D_HGW_1_l</mesh> - <field_name>p_bc</field_name> + <parameter>bc_left</parameter> </boundary_condition> <boundary_condition> - <type>NonuniformDirichlet</type> + <type>Dirichlet</type> <mesh>3D_HGW_1_r</mesh> - <field_name>p_bc</field_name> + <parameter>bc_right</parameter> </boundary_condition> </boundary_conditions> </process_variable> diff --git a/Tests/Data/ThermoMechanics/CreepBGRa/Verification/m2_2Dloadbt/m2_2Dloadbt.prj b/Tests/Data/ThermoMechanics/CreepBGRa/Verification/m2_2Dloadbt/m2_2Dloadbt.prj index 27e6952bce91e0f2fe30045fbe0eb790842845c5..61201fc6a0f0b65a2ccfb94404679bd5365bf31c 100644 --- a/Tests/Data/ThermoMechanics/CreepBGRa/Verification/m2_2Dloadbt/m2_2Dloadbt.prj +++ b/Tests/Data/ThermoMechanics/CreepBGRa/Verification/m2_2Dloadbt/m2_2Dloadbt.prj @@ -285,7 +285,7 @@ <vtkdiff> <file>m2_2Dloadbt_pcs_0_ts_300_t_0.600000.vtu</file> <field>sigma</field> - <absolute_tolerance>1e-12</absolute_tolerance> + <absolute_tolerance>5e-12</absolute_tolerance> <relative_tolerance>0</relative_tolerance> </vtkdiff> <vtkdiff>