From 7bb833c47245ec31fd84e71b399edce3ce14a5c8 Mon Sep 17 00:00:00 2001 From: Dmitri Naumov <github@naumov.de> Date: Sat, 14 May 2016 17:52:39 +0000 Subject: [PATCH] [PL] Tabs to whitespaces. --- ProcessLib/CMakeLists.txt | 4 +- ProcessLib/DirichletBc.h | 4 +- ProcessLib/InitialCondition.cpp | 62 ++-- ProcessLib/InitialCondition.h | 52 +-- ProcessLib/Parameter.cpp | 50 +-- ProcessLib/Parameter.h | 42 +-- ProcessLib/Process.cpp | 62 ++-- ProcessLib/Process.h | 594 ++++++++++++++++---------------- ProcessLib/ProcessVariable.cpp | 172 ++++----- ProcessLib/ProcessVariable.h | 118 +++---- 10 files changed, 580 insertions(+), 580 deletions(-) diff --git a/ProcessLib/CMakeLists.txt b/ProcessLib/CMakeLists.txt index a4006e51595..3a84b340bf7 100644 --- a/ProcessLib/CMakeLists.txt +++ b/ProcessLib/CMakeLists.txt @@ -18,8 +18,8 @@ target_link_libraries(ProcessLib ADD_VTK_DEPENDENCY(ProcessLib) if(TARGET Eigen) - add_dependencies(ProcessLib Eigen) + add_dependencies(ProcessLib Eigen) endif() if(TARGET Boost) - add_dependencies(ProcessLib Boost) + add_dependencies(ProcessLib Boost) endif() diff --git a/ProcessLib/DirichletBc.h b/ProcessLib/DirichletBc.h index d882b9c7916..1e608dbaa9e 100644 --- a/ProcessLib/DirichletBc.h +++ b/ProcessLib/DirichletBc.h @@ -22,8 +22,8 @@ namespace ProcessLib template <typename IndexType> struct DirichletBc final { - std::vector<IndexType> global_ids; - std::vector<double> values; + std::vector<IndexType> global_ids; + std::vector<double> values; }; } // namespace ProcessLib diff --git a/ProcessLib/InitialCondition.cpp b/ProcessLib/InitialCondition.cpp index 9058b66994f..3ea7de86908 100644 --- a/ProcessLib/InitialCondition.cpp +++ b/ProcessLib/InitialCondition.cpp @@ -23,13 +23,13 @@ namespace ProcessLib std::unique_ptr<InitialCondition> createUniformInitialCondition( BaseLib::ConfigTree const& config, int const /*n_components*/) { - config.checkConfParam("type", "Uniform"); + config.checkConfParam("type", "Uniform"); - auto value = config.getConfParam<double>("value"); - DBUG("Using value %g", value); + auto value = config.getConfParam<double>("value"); + DBUG("Using value %g", value); - return std::unique_ptr<InitialCondition>( - new UniformInitialCondition(value)); + return std::unique_ptr<InitialCondition>( + new UniformInitialCondition(value)); } std::unique_ptr<InitialCondition> createMeshPropertyInitialCondition( @@ -37,34 +37,34 @@ std::unique_ptr<InitialCondition> createMeshPropertyInitialCondition( MeshLib::Mesh const& mesh, int const n_components) { - auto field_name = config.getConfParam<std::string>("field_name"); - DBUG("Using field_name %s", field_name.c_str()); + auto field_name = config.getConfParam<std::string>("field_name"); + DBUG("Using field_name %s", field_name.c_str()); - if (!mesh.getProperties().hasPropertyVector(field_name)) - { - ERR("The required property %s does not exists in the mesh.", - field_name.c_str()); - std::abort(); - } - auto const& property = - mesh.getProperties().template getPropertyVector<double>(field_name); - if (!property) - { - ERR("The required property %s is not of the requested type.", - field_name.c_str()); - std::abort(); - } + if (!mesh.getProperties().hasPropertyVector(field_name)) + { + ERR("The required property %s does not exists in the mesh.", + field_name.c_str()); + std::abort(); + } + auto const& property = + mesh.getProperties().template getPropertyVector<double>(field_name); + if (!property) + { + ERR("The required property %s is not of the requested type.", + field_name.c_str()); + std::abort(); + } - if (property->getNumberOfComponents() != - static_cast<std::size_t>(n_components)) - { - ERR("The required property %s has different number of components %d, " - "expected %d.", - field_name.c_str(), property->getNumberOfComponents(), n_components); - std::abort(); - } - return std::unique_ptr<InitialCondition>( - new MeshPropertyInitialCondition(*property)); + if (property->getNumberOfComponents() != + static_cast<std::size_t>(n_components)) + { + ERR("The required property %s has different number of components %d, " + "expected %d.", + field_name.c_str(), property->getNumberOfComponents(), n_components); + std::abort(); + } + return std::unique_ptr<InitialCondition>( + new MeshPropertyInitialCondition(*property)); } } // namespace ProcessLib diff --git a/ProcessLib/InitialCondition.h b/ProcessLib/InitialCondition.h index a6e5400799d..bab59d97783 100644 --- a/ProcessLib/InitialCondition.h +++ b/ProcessLib/InitialCondition.h @@ -34,27 +34,27 @@ namespace ProcessLib class InitialCondition { public: - virtual ~InitialCondition() = default; - virtual double getValue(MeshLib::Node const&, int const component_id) const = 0; + virtual ~InitialCondition() = default; + virtual double getValue(MeshLib::Node const&, int const component_id) const = 0; }; /// Uniform value initial condition class UniformInitialCondition : public InitialCondition { public: - UniformInitialCondition(double const value) : _value(value) - { - } - /// Returns a value for given node and component. - /// \todo The component_id is to be implemented. - virtual double getValue(MeshLib::Node const&, - int const /* component_id */) const override - { - return _value; - } + UniformInitialCondition(double const value) : _value(value) + { + } + /// Returns a value for given node and component. + /// \todo The component_id is to be implemented. + virtual double getValue(MeshLib::Node const&, + int const /* component_id */) const override + { + return _value; + } private: - double _value; + double _value; }; /// Construct a UniformInitialCondition from configuration. @@ -68,22 +68,22 @@ std::unique_ptr<InitialCondition> createUniformInitialCondition( class MeshPropertyInitialCondition : public InitialCondition { public: - MeshPropertyInitialCondition( - MeshLib::PropertyVector<double> const& property) - : _property(property) - { - assert(_property.getMeshItemType() == MeshLib::MeshItemType::Node); - } + MeshPropertyInitialCondition( + MeshLib::PropertyVector<double> const& property) + : _property(property) + { + assert(_property.getMeshItemType() == MeshLib::MeshItemType::Node); + } - virtual double getValue(MeshLib::Node const& n, - int const component_id) const override - { - return _property[n.getID() * _property.getNumberOfComponents() + - component_id]; - } + virtual double getValue(MeshLib::Node const& n, + int const component_id) const override + { + return _property[n.getID() * _property.getNumberOfComponents() + + component_id]; + } private: - MeshLib::PropertyVector<double> const& _property; + MeshLib::PropertyVector<double> const& _property; }; /// Construct a MeshPropertyInitialCondition from configuration. diff --git a/ProcessLib/Parameter.cpp b/ProcessLib/Parameter.cpp index ecbabfe84d9..45e956acd41 100644 --- a/ProcessLib/Parameter.cpp +++ b/ProcessLib/Parameter.cpp @@ -19,36 +19,36 @@ namespace ProcessLib std::unique_ptr<ParameterBase> createConstParameter( BaseLib::ConfigTree const& config) { - config.checkConfParam("type", "Constant"); - auto value = config.getConfParam<double>("value"); - DBUG("Using value %g", value); + config.checkConfParam("type", "Constant"); + auto value = config.getConfParam<double>("value"); + DBUG("Using value %g", value); - return std::unique_ptr<ParameterBase>(new ConstParameter<double>(value)); + return std::unique_ptr<ParameterBase>(new ConstParameter<double>(value)); } std::unique_ptr<ParameterBase> createMeshPropertyParameter( BaseLib::ConfigTree const& config, MeshLib::Mesh const& mesh) { - config.checkConfParam("type", "MeshProperty"); - auto field_name = config.getConfParam<std::string>("field_name"); - DBUG("Using field_name %s", field_name.c_str()); - - if (!mesh.getProperties().hasPropertyVector(field_name)) - { - ERR("The required property %s does not exists in the mesh.", - field_name.c_str()); - std::abort(); - } - auto const& property = - mesh.getProperties().template getPropertyVector<double>(field_name); - if (!property) - { - ERR("The required property %s is not of the requested type.", - field_name.c_str()); - std::abort(); - } - - return std::unique_ptr<ParameterBase>( - new MeshPropertyParameter<double>(*property)); + config.checkConfParam("type", "MeshProperty"); + auto field_name = config.getConfParam<std::string>("field_name"); + DBUG("Using field_name %s", field_name.c_str()); + + if (!mesh.getProperties().hasPropertyVector(field_name)) + { + ERR("The required property %s does not exists in the mesh.", + field_name.c_str()); + std::abort(); + } + auto const& property = + mesh.getProperties().template getPropertyVector<double>(field_name); + if (!property) + { + ERR("The required property %s is not of the requested type.", + field_name.c_str()); + std::abort(); + } + + return std::unique_ptr<ParameterBase>( + new MeshPropertyParameter<double>(*property)); } } // namespace ProcessLib diff --git a/ProcessLib/Parameter.h b/ProcessLib/Parameter.h index 4a3d09753cd..4c52a459ee6 100644 --- a/ProcessLib/Parameter.h +++ b/ProcessLib/Parameter.h @@ -31,9 +31,9 @@ namespace ProcessLib /// Its property name helps addressing the right parameter. struct ParameterBase { - virtual ~ParameterBase() = default; + virtual ~ParameterBase() = default; - std::string name; + std::string name; }; /// A parameter is representing a value or function of any type. @@ -43,9 +43,9 @@ struct ParameterBase template <typename ReturnType, typename... Args> struct Parameter : public ParameterBase { - virtual ~Parameter() = default; + virtual ~Parameter() = default; - virtual ReturnType operator()(Args&&... args) const = 0; + virtual ReturnType operator()(Args&&... args) const = 0; }; /// Single, constant value parameter. @@ -53,17 +53,17 @@ template <typename ReturnType> struct ConstParameter final : public Parameter<ReturnType, MeshLib::Element const&> { - ConstParameter(ReturnType value) : _value(value) - { - } + ConstParameter(ReturnType value) : _value(value) + { + } - ReturnType operator()(MeshLib::Element const&) const override - { - return _value; - } + ReturnType operator()(MeshLib::Element const&) const override + { + return _value; + } private: - ReturnType _value; + ReturnType _value; }; std::unique_ptr<ParameterBase> createConstParameter(BaseLib::ConfigTree const& config); @@ -73,18 +73,18 @@ template <typename ReturnType> struct MeshPropertyParameter final : public Parameter<ReturnType, MeshLib::Element const&> { - MeshPropertyParameter(MeshLib::PropertyVector<ReturnType> const& property) - : _property(property) - { - } + MeshPropertyParameter(MeshLib::PropertyVector<ReturnType> const& property) + : _property(property) + { + } - ReturnType operator()(MeshLib::Element const& e) const override - { - return _property[e.getID()]; - } + ReturnType operator()(MeshLib::Element const& e) const override + { + return _property[e.getID()]; + } private: - MeshLib::PropertyVector<ReturnType> const& _property; + MeshLib::PropertyVector<ReturnType> const& _property; }; std::unique_ptr<ParameterBase> createMeshPropertyParameter(BaseLib::ConfigTree const& config, MeshLib::Mesh const& mesh); diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp index 556e86458b3..4ae6e359f6c 100644 --- a/ProcessLib/Process.cpp +++ b/ProcessLib/Process.cpp @@ -15,47 +15,47 @@ ProcessVariable& findProcessVariable( std::vector<ProcessVariable> const& variables, BaseLib::ConfigTree const& pv_config, std::string const& tag) { - // Find process variable name in process config. - std::string const name = pv_config.getConfParam<std::string>(tag); + // Find process variable name in process config. + std::string const name = pv_config.getConfParam<std::string>(tag); // Find corresponding variable by name. - auto variable = std::find_if(variables.cbegin(), variables.cend(), - [&name](ProcessVariable const& v) - { - return v.getName() == name; - }); - - if (variable == variables.end()) - { - ERR( - "Could not find process variable '%s' in the provided variables " - "list for config tag <%s>.", - name.c_str(), tag.c_str()); - std::abort(); - } - DBUG("Found process variable \'%s\' for config tag <%s>.", - variable->getName().c_str(), tag.c_str()); - - // Const cast is needed because of variables argument constness. - return const_cast<ProcessVariable&>(*variable); + auto variable = std::find_if(variables.cbegin(), variables.cend(), + [&name](ProcessVariable const& v) + { + return v.getName() == name; + }); + + if (variable == variables.end()) + { + ERR( + "Could not find process variable '%s' in the provided variables " + "list for config tag <%s>.", + name.c_str(), tag.c_str()); + std::abort(); + } + DBUG("Found process variable \'%s\' for config tag <%s>.", + variable->getName().c_str(), tag.c_str()); + + // Const cast is needed because of variables argument constness. + return const_cast<ProcessVariable&>(*variable); } std::vector<std::reference_wrapper<ProcessVariable>> findProcessVariables( - std::vector<ProcessVariable> const& variables, - BaseLib::ConfigTree const& process_config, - std::initializer_list<std::string> tag_names) + std::vector<ProcessVariable> const& variables, + BaseLib::ConfigTree const& process_config, + std::initializer_list<std::string> tag_names) { - std::vector<std::reference_wrapper<ProcessVariable>> vars; - vars.reserve(tag_names.size()); + std::vector<std::reference_wrapper<ProcessVariable>> vars; + vars.reserve(tag_names.size()); - auto const pv_conf = process_config.getConfSubtree("process_variables"); + auto const pv_conf = process_config.getConfSubtree("process_variables"); - for (auto const& tag : tag_names) { - vars.emplace_back(findProcessVariable(variables, pv_conf, tag)); - } + for (auto const& tag : tag_names) { + vars.emplace_back(findProcessVariable(variables, pv_conf, tag)); + } - return vars; + return vars; } } // namespace ProcessLib diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h index bda093cc0c3..c7629c60f19 100644 --- a/ProcessLib/Process.h +++ b/ProcessLib/Process.h @@ -34,294 +34,294 @@ namespace ProcessLib template <typename GlobalSetup> class Process - : public NumLib::ODESystem<typename GlobalSetup::MatrixType, - typename GlobalSetup::VectorType, - // TODO: later on use a simpler ODE system - NumLib::ODESystemTag::FirstOrderImplicitQuasilinear, - NumLib::NonlinearSolverTag::Newton> + : public NumLib::ODESystem<typename GlobalSetup::MatrixType, + typename GlobalSetup::VectorType, + // TODO: later on use a simpler ODE system + NumLib::ODESystemTag::FirstOrderImplicitQuasilinear, + NumLib::NonlinearSolverTag::Newton> { public: - using GlobalVector = typename GlobalSetup::VectorType; - using GlobalMatrix = typename GlobalSetup::MatrixType; - using Index = typename GlobalMatrix::IndexType; - using NonlinearSolver = NumLib::NonlinearSolverBase<GlobalMatrix, GlobalVector>; - using TimeDiscretization = NumLib::TimeDiscretization<GlobalVector>; - - Process( - MeshLib::Mesh& mesh, - NonlinearSolver& nonlinear_solver, - std::unique_ptr<TimeDiscretization>&& time_discretization, - std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables, - SecondaryVariableCollection<GlobalVector>&& secondary_variables, - ProcessOutput<GlobalVector>&& process_output - ) - : _mesh(mesh) - , _nonlinear_solver(nonlinear_solver) - , _time_discretization(std::move(time_discretization)) - , _process_variables(std::move(process_variables)) - , _process_output(std::move(process_output)) - , _secondary_variables(std::move(secondary_variables)) - {} - - /// Preprocessing before starting assembly for new timestep. - virtual void preTimestep(GlobalVector const& /*x*/, - const double /*t*/, const double /*delta_t*/) {} - - /// Postprocessing after a complete timestep. - virtual void postTimestep(GlobalVector const& /*x*/) {} - - /// Process output. - /// The file_name is indicating the name of possible output file. - void output(std::string const& file_name, - const unsigned /*timestep*/, - GlobalVector const& x) const - { - doProcessOutput(file_name, x, _mesh, *_local_to_global_index_map, - _process_variables, _secondary_variables, _process_output); - } - - void initialize() - { - DBUG("Initialize process."); - - DBUG("Construct dof mappings."); - constructDofTable(); + using GlobalVector = typename GlobalSetup::VectorType; + using GlobalMatrix = typename GlobalSetup::MatrixType; + using Index = typename GlobalMatrix::IndexType; + using NonlinearSolver = NumLib::NonlinearSolverBase<GlobalMatrix, GlobalVector>; + using TimeDiscretization = NumLib::TimeDiscretization<GlobalVector>; + + Process( + MeshLib::Mesh& mesh, + NonlinearSolver& nonlinear_solver, + std::unique_ptr<TimeDiscretization>&& time_discretization, + std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables, + SecondaryVariableCollection<GlobalVector>&& secondary_variables, + ProcessOutput<GlobalVector>&& process_output + ) + : _mesh(mesh) + , _nonlinear_solver(nonlinear_solver) + , _time_discretization(std::move(time_discretization)) + , _process_variables(std::move(process_variables)) + , _process_output(std::move(process_output)) + , _secondary_variables(std::move(secondary_variables)) + {} + + /// Preprocessing before starting assembly for new timestep. + virtual void preTimestep(GlobalVector const& /*x*/, + const double /*t*/, const double /*delta_t*/) {} + + /// Postprocessing after a complete timestep. + virtual void postTimestep(GlobalVector const& /*x*/) {} + + /// Process output. + /// The file_name is indicating the name of possible output file. + void output(std::string const& file_name, + const unsigned /*timestep*/, + GlobalVector const& x) const + { + doProcessOutput(file_name, x, _mesh, *_local_to_global_index_map, + _process_variables, _secondary_variables, _process_output); + } + + void initialize() + { + DBUG("Initialize process."); + + DBUG("Construct dof mappings."); + constructDofTable(); #ifndef USE_PETSC - DBUG("Compute sparsity pattern"); - computeSparsityPattern(); + DBUG("Compute sparsity pattern"); + computeSparsityPattern(); #endif - initializeConcreteProcess(*_local_to_global_index_map, _mesh, - _integration_order); - - DBUG("Initialize boundary conditions."); - for (ProcessVariable& pv : _process_variables) - { - createDirichletBcs(pv, 0); // 0 is the component id - createNeumannBcs(pv, 0); // 0 is the component id - } - - for (auto& bc : _neumann_bcs) - bc->initialize(_mesh.getDimension()); - } - - void setInitialConditions(GlobalVector& x) - { - DBUG("Set initial conditions."); - for (ProcessVariable& pv : _process_variables) - { - setInitialConditions(pv, 0, x); // 0 is the component id - } - } - - MathLib::MatrixSpecifications getMatrixSpecifications() const override final - { - return { 0u, 0u, - &_sparsity_pattern, - _local_to_global_index_map.get(), - &_mesh }; - } - - void assemble(const double t, GlobalVector const& x, - GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) override final - { - assembleConcreteProcess(t, x, M, K, b); - - // Call global assembler for each Neumann boundary local assembler. - for (auto const& bc : _neumann_bcs) - bc->integrate(t, b); - } - - void assembleJacobian( - const double t, GlobalVector const& x, GlobalVector const& xdot, - const double dxdot_dx, GlobalMatrix const& M, - const double dx_dx, GlobalMatrix const& K, - GlobalMatrix& Jac) override final - { - assembleJacobianConcreteProcess(t, x, xdot, dxdot_dx, M, dx_dx, K, Jac); - - // TODO In this method one could check if the user wants to use an - // analytical or a numerical Jacobian. Then the right - // assembleJacobianConcreteProcess() method will be chosen. - // Additionally in the default implementation of said method one - // could provide a fallback to a numerical Jacobian. However, that - // would be in a sense implicit behaviour and it might be better to - // just abort, as is currently the case. - // In order to implement the Jacobian assembly entirely, in addition - // to the operator() in VectorMatrixAssembler there has to be a method - // that dispatches the Jacobian assembly. - // Similarly, then the NeumannBC specialization of VectorMatrixAssembler - // probably can be merged into the main class s.t. one has only one - // type of VectorMatrixAssembler (for each equation type) with the - // three methods assemble(), assembleJacobian() and assembleNeumannBC(). - // That list can be extended, e.g. by methods for the assembly of - // source terms. - // UPDATE: Probably it is better to keep a separate NeumannBC version of the - // VectoMatrixAssembler since that will work for all kinds of processes. - } - - std::vector<DirichletBc<Index>> const* getKnownSolutions( - double const /*t*/) const override final - { - return &_dirichlet_bcs; - } - - NonlinearSolver& getNonlinearSolver() const - { - return _nonlinear_solver; - } - - TimeDiscretization& getTimeDiscretization() const - { - return *_time_discretization; - } + initializeConcreteProcess(*_local_to_global_index_map, _mesh, + _integration_order); + + DBUG("Initialize boundary conditions."); + for (ProcessVariable& pv : _process_variables) + { + createDirichletBcs(pv, 0); // 0 is the component id + createNeumannBcs(pv, 0); // 0 is the component id + } + + for (auto& bc : _neumann_bcs) + bc->initialize(_mesh.getDimension()); + } + + void setInitialConditions(GlobalVector& x) + { + DBUG("Set initial conditions."); + for (ProcessVariable& pv : _process_variables) + { + setInitialConditions(pv, 0, x); // 0 is the component id + } + } + + MathLib::MatrixSpecifications getMatrixSpecifications() const override final + { + return { 0u, 0u, + &_sparsity_pattern, + _local_to_global_index_map.get(), + &_mesh }; + } + + void assemble(const double t, GlobalVector const& x, + GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) override final + { + assembleConcreteProcess(t, x, M, K, b); + + // Call global assembler for each Neumann boundary local assembler. + for (auto const& bc : _neumann_bcs) + bc->integrate(t, b); + } + + void assembleJacobian( + const double t, GlobalVector const& x, GlobalVector const& xdot, + const double dxdot_dx, GlobalMatrix const& M, + const double dx_dx, GlobalMatrix const& K, + GlobalMatrix& Jac) override final + { + assembleJacobianConcreteProcess(t, x, xdot, dxdot_dx, M, dx_dx, K, Jac); + + // TODO In this method one could check if the user wants to use an + // analytical or a numerical Jacobian. Then the right + // assembleJacobianConcreteProcess() method will be chosen. + // Additionally in the default implementation of said method one + // could provide a fallback to a numerical Jacobian. However, that + // would be in a sense implicit behaviour and it might be better to + // just abort, as is currently the case. + // In order to implement the Jacobian assembly entirely, in addition + // to the operator() in VectorMatrixAssembler there has to be a method + // that dispatches the Jacobian assembly. + // Similarly, then the NeumannBC specialization of VectorMatrixAssembler + // probably can be merged into the main class s.t. one has only one + // type of VectorMatrixAssembler (for each equation type) with the + // three methods assemble(), assembleJacobian() and assembleNeumannBC(). + // That list can be extended, e.g. by methods for the assembly of + // source terms. + // UPDATE: Probably it is better to keep a separate NeumannBC version of the + // VectoMatrixAssembler since that will work for all kinds of processes. + } + + std::vector<DirichletBc<Index>> const* getKnownSolutions( + double const /*t*/) const override final + { + return &_dirichlet_bcs; + } + + NonlinearSolver& getNonlinearSolver() const + { + return _nonlinear_solver; + } + + TimeDiscretization& getTimeDiscretization() const + { + return *_time_discretization; + } private: - /// Process specific initialization called by initialize(). - virtual void initializeConcreteProcess( - AssemblerLib::LocalToGlobalIndexMap const& dof_table, - MeshLib::Mesh const& mesh, - unsigned const integration_order) = 0; - - virtual void assembleConcreteProcess( - const double t, GlobalVector const& x, - GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) = 0; - - virtual void assembleJacobianConcreteProcess( - const double /*t*/, GlobalVector const& /*x*/, GlobalVector const& /*xdot*/, - const double /*dxdot_dx*/, GlobalMatrix const& /*M*/, - const double /*dx_dx*/, GlobalMatrix const& /*K*/, - GlobalMatrix& /*Jac*/) - { - ERR("The concrete implementation of this Process did not override the" - " assembleJacobianConcreteProcess() method." - " Hence, no analytical Jacobian is provided for this process" - " and the Newton-Raphson method cannot be used to solve it."); - std::abort(); - } - - void constructDofTable() - { - // Create single component dof in every of the mesh's nodes. - _mesh_subset_all_nodes.reset( - new MeshLib::MeshSubset(_mesh, &_mesh.getNodes())); - - // Collect the mesh subsets in a vector. - std::vector<std::unique_ptr<MeshLib::MeshSubsets>> all_mesh_subsets; - for (ProcessVariable const& pv : _process_variables) - { - std::generate_n( - std::back_inserter(all_mesh_subsets), - pv.getNumberOfComponents(), - [&]() - { - return std::unique_ptr<MeshLib::MeshSubsets>{ - new MeshLib::MeshSubsets{_mesh_subset_all_nodes.get()}}; - }); - } - - _local_to_global_index_map.reset( - new AssemblerLib::LocalToGlobalIndexMap( - std::move(all_mesh_subsets), - AssemblerLib::ComponentOrder::BY_LOCATION)); - } - - /// Sets the initial condition values in the solution vector x for a given - /// process variable and component. - void setInitialConditions(ProcessVariable const& variable, - int const component_id, - GlobalVector& x) - { - std::size_t const n_nodes = _mesh.getNNodes(); - for (std::size_t node_id = 0; node_id < n_nodes; ++node_id) - { - MeshLib::Location const l(_mesh.getID(), - MeshLib::MeshItemType::Node, node_id); - auto global_index = std::abs( - _local_to_global_index_map->getGlobalIndex(l, component_id)); + /// Process specific initialization called by initialize(). + virtual void initializeConcreteProcess( + AssemblerLib::LocalToGlobalIndexMap const& dof_table, + MeshLib::Mesh const& mesh, + unsigned const integration_order) = 0; + + virtual void assembleConcreteProcess( + const double t, GlobalVector const& x, + GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) = 0; + + virtual void assembleJacobianConcreteProcess( + const double /*t*/, GlobalVector const& /*x*/, GlobalVector const& /*xdot*/, + const double /*dxdot_dx*/, GlobalMatrix const& /*M*/, + const double /*dx_dx*/, GlobalMatrix const& /*K*/, + GlobalMatrix& /*Jac*/) + { + ERR("The concrete implementation of this Process did not override the" + " assembleJacobianConcreteProcess() method." + " Hence, no analytical Jacobian is provided for this process" + " and the Newton-Raphson method cannot be used to solve it."); + std::abort(); + } + + void constructDofTable() + { + // Create single component dof in every of the mesh's nodes. + _mesh_subset_all_nodes.reset( + new MeshLib::MeshSubset(_mesh, &_mesh.getNodes())); + + // Collect the mesh subsets in a vector. + std::vector<std::unique_ptr<MeshLib::MeshSubsets>> all_mesh_subsets; + for (ProcessVariable const& pv : _process_variables) + { + std::generate_n( + std::back_inserter(all_mesh_subsets), + pv.getNumberOfComponents(), + [&]() + { + return std::unique_ptr<MeshLib::MeshSubsets>{ + new MeshLib::MeshSubsets{_mesh_subset_all_nodes.get()}}; + }); + } + + _local_to_global_index_map.reset( + new AssemblerLib::LocalToGlobalIndexMap( + std::move(all_mesh_subsets), + AssemblerLib::ComponentOrder::BY_LOCATION)); + } + + /// Sets the initial condition values in the solution vector x for a given + /// process variable and component. + void setInitialConditions(ProcessVariable const& variable, + int const component_id, + GlobalVector& x) + { + std::size_t const n_nodes = _mesh.getNNodes(); + for (std::size_t node_id = 0; node_id < n_nodes; ++node_id) + { + MeshLib::Location const l(_mesh.getID(), + MeshLib::MeshItemType::Node, node_id); + auto global_index = std::abs( + _local_to_global_index_map->getGlobalIndex(l, component_id)); #ifdef USE_PETSC - // The global indices of the ghost entries of the global - // matrix or the global vectors need to be set as negative values - // for equation assembly, however the global indices start from zero. - // Therefore, any ghost entry with zero index is assigned an negative - // value of the vector size or the matrix dimension. - // To assign the initial value for the ghost entries, - // the negative indices of the ghost entries are restored to zero. - // checked hereby. - if ( global_index == x.size() ) - global_index = 0; + // The global indices of the ghost entries of the global + // matrix or the global vectors need to be set as negative values + // for equation assembly, however the global indices start from zero. + // Therefore, any ghost entry with zero index is assigned an negative + // value of the vector size or the matrix dimension. + // To assign the initial value for the ghost entries, + // the negative indices of the ghost entries are restored to zero. + // checked hereby. + if ( global_index == x.size() ) + global_index = 0; #endif - x.set(global_index, - variable.getInitialConditionValue(*_mesh.getNode(node_id), - component_id)); - } - } - - void createDirichletBcs(ProcessVariable& variable, int const component_id) - { - MeshGeoToolsLib::MeshNodeSearcher& mesh_node_searcher = - MeshGeoToolsLib::MeshNodeSearcher::getMeshNodeSearcher( - variable.getMesh()); - - variable.initializeDirichletBCs(std::back_inserter(_dirichlet_bcs), - mesh_node_searcher, - *_local_to_global_index_map, - component_id); - } - - void createNeumannBcs(ProcessVariable& variable, int const component_id) - { - // Find mesh nodes. - MeshGeoToolsLib::MeshNodeSearcher& mesh_node_searcher = - MeshGeoToolsLib::MeshNodeSearcher::getMeshNodeSearcher( - variable.getMesh()); - MeshGeoToolsLib::BoundaryElementsSearcher mesh_element_searcher( - variable.getMesh(), mesh_node_searcher); - - // Create a neumann BC for the process variable storing them in the - // _neumann_bcs vector. - variable.createNeumannBcs<GlobalSetup>( - std::back_inserter(_neumann_bcs), - mesh_element_searcher, - _integration_order, - *_local_to_global_index_map, - 0, // 0 is the variable id TODO - component_id); - } - - /// Computes and stores global matrix' sparsity pattern from given - /// DOF-table. - void computeSparsityPattern() - { - _sparsity_pattern = std::move(AssemblerLib::computeSparsityPattern( - *_local_to_global_index_map, _mesh)); - } + x.set(global_index, + variable.getInitialConditionValue(*_mesh.getNode(node_id), + component_id)); + } + } + + void createDirichletBcs(ProcessVariable& variable, int const component_id) + { + MeshGeoToolsLib::MeshNodeSearcher& mesh_node_searcher = + MeshGeoToolsLib::MeshNodeSearcher::getMeshNodeSearcher( + variable.getMesh()); + + variable.initializeDirichletBCs(std::back_inserter(_dirichlet_bcs), + mesh_node_searcher, + *_local_to_global_index_map, + component_id); + } + + void createNeumannBcs(ProcessVariable& variable, int const component_id) + { + // Find mesh nodes. + MeshGeoToolsLib::MeshNodeSearcher& mesh_node_searcher = + MeshGeoToolsLib::MeshNodeSearcher::getMeshNodeSearcher( + variable.getMesh()); + MeshGeoToolsLib::BoundaryElementsSearcher mesh_element_searcher( + variable.getMesh(), mesh_node_searcher); + + // Create a neumann BC for the process variable storing them in the + // _neumann_bcs vector. + variable.createNeumannBcs<GlobalSetup>( + std::back_inserter(_neumann_bcs), + mesh_element_searcher, + _integration_order, + *_local_to_global_index_map, + 0, // 0 is the variable id TODO + component_id); + } + + /// Computes and stores global matrix' sparsity pattern from given + /// DOF-table. + void computeSparsityPattern() + { + _sparsity_pattern = std::move(AssemblerLib::computeSparsityPattern( + *_local_to_global_index_map, _mesh)); + } private: - unsigned const _integration_order = 2; + unsigned const _integration_order = 2; - MeshLib::Mesh& _mesh; + MeshLib::Mesh& _mesh; - std::unique_ptr<AssemblerLib::LocalToGlobalIndexMap> - _local_to_global_index_map; + std::unique_ptr<AssemblerLib::LocalToGlobalIndexMap> + _local_to_global_index_map; - MathLib::SparsityPattern _sparsity_pattern; + MathLib::SparsityPattern _sparsity_pattern; - std::vector<DirichletBc<GlobalIndexType>> _dirichlet_bcs; - std::vector<std::unique_ptr<NeumannBc<GlobalSetup>>> _neumann_bcs; + std::vector<DirichletBc<GlobalIndexType>> _dirichlet_bcs; + std::vector<std::unique_ptr<NeumannBc<GlobalSetup>>> _neumann_bcs; - NonlinearSolver& _nonlinear_solver; - std::unique_ptr<TimeDiscretization> _time_discretization; + NonlinearSolver& _nonlinear_solver; + std::unique_ptr<TimeDiscretization> _time_discretization; - /// Variables used by this process. - std::vector<std::reference_wrapper<ProcessVariable>> _process_variables; + /// Variables used by this process. + std::vector<std::reference_wrapper<ProcessVariable>> _process_variables; - ProcessOutput<GlobalVector> _process_output; + ProcessOutput<GlobalVector> _process_output; protected: - std::unique_ptr<MeshLib::MeshSubset const> _mesh_subset_all_nodes; - SecondaryVariableCollection<GlobalVector> _secondary_variables; + std::unique_ptr<MeshLib::MeshSubset const> _mesh_subset_all_nodes; + SecondaryVariableCollection<GlobalVector> _secondary_variables; }; /// Find process variables in \c variables whose names match the settings under @@ -366,36 +366,36 @@ Parameter<ParameterArgs...>& findParameter( BaseLib::ConfigTree const& process_config, std::string const& tag, std::vector<std::unique_ptr<ParameterBase>> const& parameters) { - // Find parameter name in process config. - auto const name = process_config.getConfParam<std::string>(tag); - - // Find corresponding parameter by name. - auto const parameter_it = - std::find_if(parameters.cbegin(), parameters.cend(), - [&name](std::unique_ptr<ParameterBase> const& p) - { - return p->name == name; - }); - - if (parameter_it == parameters.end()) - { - ERR( - "Could not find parameter '%s' in the provided parameters list for " - "config tag <%s>.", - name.c_str(), tag.c_str()); - std::abort(); - } - DBUG("Found parameter \'%s\'.", (*parameter_it)->name.c_str()); - - // Check the type correctness of the found parameter. - auto* const parameter = - dynamic_cast<Parameter<ParameterArgs...>*>(parameter_it->get()); - if (!parameter) - { - ERR("The read parameter is of incompatible type."); - std::abort(); - } - return *parameter; + // Find parameter name in process config. + auto const name = process_config.getConfParam<std::string>(tag); + + // Find corresponding parameter by name. + auto const parameter_it = + std::find_if(parameters.cbegin(), parameters.cend(), + [&name](std::unique_ptr<ParameterBase> const& p) + { + return p->name == name; + }); + + if (parameter_it == parameters.end()) + { + ERR( + "Could not find parameter '%s' in the provided parameters list for " + "config tag <%s>.", + name.c_str(), tag.c_str()); + std::abort(); + } + DBUG("Found parameter \'%s\'.", (*parameter_it)->name.c_str()); + + // Check the type correctness of the found parameter. + auto* const parameter = + dynamic_cast<Parameter<ParameterArgs...>*>(parameter_it->get()); + if (!parameter) + { + ERR("The read parameter is of incompatible type."); + std::abort(); + } + return *parameter; } } // namespace ProcessLib diff --git a/ProcessLib/ProcessVariable.cpp b/ProcessLib/ProcessVariable.cpp index d6466c91043..953f2800453 100644 --- a/ProcessLib/ProcessVariable.cpp +++ b/ProcessLib/ProcessVariable.cpp @@ -25,74 +25,74 @@ ProcessVariable::ProcessVariable(BaseLib::ConfigTree const& config, _mesh(mesh), _n_components(config.getConfParam<int>("components")) { - DBUG("Constructing process variable %s", this->_name.c_str()); - - // Initial condition - if (auto ic_config = config.getConfSubtreeOptional("initial_condition")) - { - auto const type = ic_config->peekConfParam<std::string>("type"); - if (type == "Uniform") - { - _initial_condition = - createUniformInitialCondition(*ic_config, _n_components); - } - else if (type == "MeshProperty") - { - _initial_condition = - createMeshPropertyInitialCondition(*ic_config, _mesh, _n_components); - } - else - { - ERR("Unknown type of the initial condition."); - } - } - else - { - INFO("No initial condition found."); - } - - // Boundary conditions - if (auto bcs_config = config.getConfSubtreeOptional("boundary_conditions")) - { - for (auto bc_config - : bcs_config->getConfSubtreeList("boundary_condition")) - { - auto const geometrical_set_name = - bc_config.getConfParam<std::string>("geometrical_set"); - auto const geometry_name = - bc_config.getConfParam<std::string>("geometry"); - - GeoLib::GeoObject const* const geometry = - geometries.getGeoObject(geometrical_set_name, geometry_name); - DBUG( - "Found geometry type \"%s\"", - GeoLib::convertGeoTypeToString(geometry->getGeoType()).c_str()); - - // Construct type dependent boundary condition - auto const type = bc_config.peekConfParam<std::string>("type"); - - if (type == "UniformDirichlet") - { - _dirichlet_bc_configs.emplace_back( - new UniformDirichletBoundaryCondition(geometry, bc_config)); - } - else if (type == "UniformNeumann") - { - _neumann_bc_configs.emplace_back( - new NeumannBcConfig(geometry, bc_config)); - } - else - { - ERR("Unknown type \'%s\' of the boundary condition.", - type.c_str()); - } - } - } else { - INFO("No boundary conditions found."); - } - - // Source Terms - config.ignoreConfParam("source_terms"); + DBUG("Constructing process variable %s", this->_name.c_str()); + + // Initial condition + if (auto ic_config = config.getConfSubtreeOptional("initial_condition")) + { + auto const type = ic_config->peekConfParam<std::string>("type"); + if (type == "Uniform") + { + _initial_condition = + createUniformInitialCondition(*ic_config, _n_components); + } + else if (type == "MeshProperty") + { + _initial_condition = + createMeshPropertyInitialCondition(*ic_config, _mesh, _n_components); + } + else + { + ERR("Unknown type of the initial condition."); + } + } + else + { + INFO("No initial condition found."); + } + + // Boundary conditions + if (auto bcs_config = config.getConfSubtreeOptional("boundary_conditions")) + { + for (auto bc_config + : bcs_config->getConfSubtreeList("boundary_condition")) + { + auto const geometrical_set_name = + bc_config.getConfParam<std::string>("geometrical_set"); + auto const geometry_name = + bc_config.getConfParam<std::string>("geometry"); + + GeoLib::GeoObject const* const geometry = + geometries.getGeoObject(geometrical_set_name, geometry_name); + DBUG( + "Found geometry type \"%s\"", + GeoLib::convertGeoTypeToString(geometry->getGeoType()).c_str()); + + // Construct type dependent boundary condition + auto const type = bc_config.peekConfParam<std::string>("type"); + + if (type == "UniformDirichlet") + { + _dirichlet_bc_configs.emplace_back( + new UniformDirichletBoundaryCondition(geometry, bc_config)); + } + else if (type == "UniformNeumann") + { + _neumann_bc_configs.emplace_back( + new NeumannBcConfig(geometry, bc_config)); + } + else + { + ERR("Unknown type \'%s\' of the boundary condition.", + type.c_str()); + } + } + } else { + INFO("No boundary conditions found."); + } + + // Source Terms + config.ignoreConfParam("source_terms"); } ProcessVariable::ProcessVariable(ProcessVariable&& other) @@ -107,32 +107,32 @@ ProcessVariable::ProcessVariable(ProcessVariable&& other) std::string const& ProcessVariable::getName() const { - return _name; + return _name; } MeshLib::Mesh const& ProcessVariable::getMesh() const { - return _mesh; + return _mesh; } MeshLib::PropertyVector<double>& ProcessVariable::getOrCreateMeshProperty() { - boost::optional<MeshLib::PropertyVector<double>&> result; - if (_mesh.getProperties().hasPropertyVector(_name)) - { - result = - _mesh.getProperties().template getPropertyVector<double>(_name); - assert(result); - assert(result->size() == _mesh.getNNodes() * _n_components); - } - else - { - result = _mesh.getProperties().template createNewPropertyVector<double>( - _name, MeshLib::MeshItemType::Node); - assert(result); - result->resize(_mesh.getNNodes() * _n_components); - } - return *result; + boost::optional<MeshLib::PropertyVector<double>&> result; + if (_mesh.getProperties().hasPropertyVector(_name)) + { + result = + _mesh.getProperties().template getPropertyVector<double>(_name); + assert(result); + assert(result->size() == _mesh.getNNodes() * _n_components); + } + else + { + result = _mesh.getProperties().template createNewPropertyVector<double>( + _name, MeshLib::MeshItemType::Node); + assert(result); + result->resize(_mesh.getNNodes() * _n_components); + } + return *result; } } // namespace ProcessLib diff --git a/ProcessLib/ProcessVariable.h b/ProcessLib/ProcessVariable.h index 66585b2b05c..6386cbb396e 100644 --- a/ProcessLib/ProcessVariable.h +++ b/ProcessLib/ProcessVariable.h @@ -47,67 +47,67 @@ namespace ProcessLib class ProcessVariable { public: - ProcessVariable(BaseLib::ConfigTree const& config, MeshLib::Mesh& mesh, - GeoLib::GEOObjects const& geometries); - - ProcessVariable(ProcessVariable&&); - - std::string const& getName() const; - - /// Returns a mesh on which the process variable is defined. - MeshLib::Mesh const& getMesh() const; - - /// Returns the number of components of the process variable. - int getNumberOfComponents() const { return _n_components; } - - template <typename OutputIterator> - void initializeDirichletBCs( - OutputIterator output_bcs, - MeshGeoToolsLib::MeshNodeSearcher& searcher, - const AssemblerLib::LocalToGlobalIndexMap& dof_table, - const unsigned nodal_dof_idx) - { - for (auto& bc_config : _dirichlet_bc_configs) - { - DirichletBc<GlobalIndexType> bc; - bc_config->initialize(searcher, dof_table, nodal_dof_idx, bc); - output_bcs++ = bc; - } - } - - template <typename GlobalSetup, typename OutputIterator, typename... Args> - void createNeumannBcs(OutputIterator bcs, - MeshGeoToolsLib::BoundaryElementsSearcher& searcher, - Args&&... args) - { - for (auto& config : _neumann_bc_configs) - { - config->initialize(searcher); - bcs++ = std::unique_ptr<NeumannBc<GlobalSetup>>{ - new NeumannBc<GlobalSetup>(*config, - std::forward<Args>(args)...)}; - } - } - - double getInitialConditionValue(MeshLib::Node const& n, - int const component_id) const - { - return _initial_condition->getValue(n, component_id); - } - - // Get or create a property vector for results. - // The returned mesh property size is number of mesh nodes times number of - // components. - MeshLib::PropertyVector<double>& getOrCreateMeshProperty(); + ProcessVariable(BaseLib::ConfigTree const& config, MeshLib::Mesh& mesh, + GeoLib::GEOObjects const& geometries); + + ProcessVariable(ProcessVariable&&); + + std::string const& getName() const; + + /// Returns a mesh on which the process variable is defined. + MeshLib::Mesh const& getMesh() const; + + /// Returns the number of components of the process variable. + int getNumberOfComponents() const { return _n_components; } + + template <typename OutputIterator> + void initializeDirichletBCs( + OutputIterator output_bcs, + MeshGeoToolsLib::MeshNodeSearcher& searcher, + const AssemblerLib::LocalToGlobalIndexMap& dof_table, + const unsigned nodal_dof_idx) + { + for (auto& bc_config : _dirichlet_bc_configs) + { + DirichletBc<GlobalIndexType> bc; + bc_config->initialize(searcher, dof_table, nodal_dof_idx, bc); + output_bcs++ = bc; + } + } + + template <typename GlobalSetup, typename OutputIterator, typename... Args> + void createNeumannBcs(OutputIterator bcs, + MeshGeoToolsLib::BoundaryElementsSearcher& searcher, + Args&&... args) + { + for (auto& config : _neumann_bc_configs) + { + config->initialize(searcher); + bcs++ = std::unique_ptr<NeumannBc<GlobalSetup>>{ + new NeumannBc<GlobalSetup>(*config, + std::forward<Args>(args)...)}; + } + } + + double getInitialConditionValue(MeshLib::Node const& n, + int const component_id) const + { + return _initial_condition->getValue(n, component_id); + } + + // Get or create a property vector for results. + // The returned mesh property size is number of mesh nodes times number of + // components. + MeshLib::PropertyVector<double>& getOrCreateMeshProperty(); private: - std::string const _name; - MeshLib::Mesh& _mesh; - const int _n_components; - std::unique_ptr<InitialCondition> _initial_condition; - std::vector<std::unique_ptr<UniformDirichletBoundaryCondition>> - _dirichlet_bc_configs; - std::vector<std::unique_ptr<NeumannBcConfig>> _neumann_bc_configs; + std::string const _name; + MeshLib::Mesh& _mesh; + const int _n_components; + std::unique_ptr<InitialCondition> _initial_condition; + std::vector<std::unique_ptr<UniformDirichletBoundaryCondition>> + _dirichlet_bc_configs; + std::vector<std::unique_ptr<NeumannBcConfig>> _neumann_bc_configs; }; } // namespace ProcessLib -- GitLab