diff --git a/ProcessLib/HeatConduction/CreateHeatConductionProcess.cpp b/ProcessLib/HeatConduction/CreateHeatConductionProcess.cpp index 8241a0b33a1f0acdbdc0e39dec8aabf237cca0ef..47fa8a0c1b08069e85b957e61d2ccbc5ea7625c9 100644 --- a/ProcessLib/HeatConduction/CreateHeatConductionProcess.cpp +++ b/ProcessLib/HeatConduction/CreateHeatConductionProcess.cpp @@ -20,6 +20,7 @@ namespace HeatConduction { std::unique_ptr<Process> createHeatConductionProcess( MeshLib::Mesh& mesh, + std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler, std::vector<ProcessVariable> const& variables, std::vector<std::unique_ptr<ParameterBase>> const& parameters, BaseLib::ConfigTree const& config) @@ -72,7 +73,8 @@ std::unique_ptr<Process> createHeatConductionProcess( named_function_caller); return std::unique_ptr<Process>{new HeatConductionProcess{ - mesh, parameters, std::move(process_variables), std::move(process_data), + mesh, std::move(jacobian_assembler), parameters, + std::move(process_variables), std::move(process_data), std::move(secondary_variables), std::move(named_function_caller)}}; } diff --git a/ProcessLib/HeatConduction/CreateHeatConductionProcess.h b/ProcessLib/HeatConduction/CreateHeatConductionProcess.h index 43979abaa7733645de5e06d4f70d0c0ef45ec0dc..36efe857499e7c79c1d0e0d30be8c71820421348 100644 --- a/ProcessLib/HeatConduction/CreateHeatConductionProcess.h +++ b/ProcessLib/HeatConduction/CreateHeatConductionProcess.h @@ -19,6 +19,7 @@ namespace HeatConduction { std::unique_ptr<Process> createHeatConductionProcess( MeshLib::Mesh& mesh, + std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler, std::vector<ProcessVariable> const& variables, std::vector<std::unique_ptr<ParameterBase>> const& parameters, BaseLib::ConfigTree const& config); diff --git a/ProcessLib/HeatConduction/HeatConductionFEM.h b/ProcessLib/HeatConduction/HeatConductionFEM.h index 5289d59e9685dcdd9ee244f9d3213037c278eadc..f4cdd26d1861710a63e1b12ac953ad1258623797 100644 --- a/ProcessLib/HeatConduction/HeatConductionFEM.h +++ b/ProcessLib/HeatConduction/HeatConductionFEM.h @@ -12,6 +12,8 @@ #include <vector> +#include "HeatConductionProcessData.h" +#include "MathLib/LinAlg/Eigen/EigenMapTools.h" #include "NumLib/Extrapolation/ExtrapolatableElement.h" #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h" #include "NumLib/Fem/ShapeMatrixPolicy.h" @@ -19,7 +21,6 @@ #include "ProcessLib/LocalAssemblerTraits.h" #include "ProcessLib/Parameter/Parameter.h" #include "ProcessLib/Utils/InitShapeMatrices.h" -#include "HeatConductionProcessData.h" namespace ProcessLib { @@ -65,10 +66,6 @@ public: HeatConductionProcessData const& process_data) : _element(element), _process_data(process_data), - _localK(local_matrix_size, - local_matrix_size), // TODO narrowing conversion - _localM(local_matrix_size, local_matrix_size), - _localRhs(local_matrix_size), _integration_method(integration_order), _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType, IntegrationMethod, GlobalDim>( @@ -79,14 +76,19 @@ public: assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF); } - void assembleConcrete( - double const t, std::vector<double> const& local_x, - NumLib::LocalToGlobalIndexMap::RowColumnIndices const& indices, - GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) override + void assemble(double const t, std::vector<double> const& local_x, + std::vector<double>& local_M_data, + std::vector<double>& local_K_data, + std::vector<double>& /*local_b_data*/) override { - _localK.setZero(); - _localM.setZero(); - _localRhs.setZero(); + auto const local_matrix_size = local_x.size(); + // This assertion is valid only if all nodal d.o.f. use the same shape matrices. + assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF); + + auto local_M = MathLib::createZeroedMatrix<NodalMatrixType>( + local_M_data, local_matrix_size, local_matrix_size); + auto local_K = MathLib::createZeroedMatrix<NodalMatrixType>( + local_K_data, local_matrix_size, local_matrix_size); unsigned const n_integration_points = _integration_method.getNumberOfPoints(); @@ -103,9 +105,9 @@ public: auto const heat_capacity = _process_data.heat_capacity(t, pos)[0]; auto const density = _process_data.density(t, pos)[0]; - _localK.noalias() += + local_K.noalias() += sm.dNdx.transpose() * k * sm.dNdx * sm.detJ * wp.getWeight(); - _localM.noalias() += sm.N.transpose() * density * heat_capacity * + local_M.noalias() += sm.N.transpose() * density * heat_capacity * sm.N * sm.detJ * wp.getWeight(); // heat flux only computed for output. GlobalDimVectorType const heat_flux = @@ -117,10 +119,6 @@ public: _heat_fluxes[d][ip] = heat_flux[d]; } } - - K.add(indices, _localK); - M.add(indices, _localM); - b.add(indices.rows, _localRhs); } Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix( @@ -157,10 +155,6 @@ private: MeshLib::Element const& _element; HeatConductionProcessData const& _process_data; - NodalMatrixType _localK; - NodalMatrixType _localM; - NodalVectorType _localRhs; - IntegrationMethod const _integration_method; std::vector<ShapeMatrices> _shape_matrices; diff --git a/ProcessLib/HeatConduction/HeatConductionProcess.cpp b/ProcessLib/HeatConduction/HeatConductionProcess.cpp index 1152cd0db1a24e562bc7a5244b09b93b607d1e2d..8523f8648701fe0ccf280066e2295343c6dd902e 100644 --- a/ProcessLib/HeatConduction/HeatConductionProcess.cpp +++ b/ProcessLib/HeatConduction/HeatConductionProcess.cpp @@ -19,13 +19,15 @@ namespace HeatConduction { HeatConductionProcess::HeatConductionProcess( MeshLib::Mesh& mesh, + std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler, std::vector<std::unique_ptr<ParameterBase>> const& parameters, std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables, HeatConductionProcessData&& process_data, SecondaryVariableCollection&& secondary_variables, NumLib::NamedFunctionCaller&& named_function_caller) - : Process(mesh, parameters, std::move(process_variables), - std::move(secondary_variables), std::move(named_function_caller)), + : Process(mesh, std::move(jacobian_assembler), parameters, + std::move(process_variables), std::move(secondary_variables), + std::move(named_function_caller)), _process_data(std::move(process_data)) { } @@ -72,10 +74,24 @@ void HeatConductionProcess::assembleConcreteProcess(const double t, DBUG("Assemble HeatConductionProcess."); // Call global assembler for each local assembly item. - GlobalExecutor::executeMemberOnDereferenced( - &HeatConductionLocalAssemblerInterface::assemble, _local_assemblers, + GlobalExecutor::executeMemberDereferenced( + _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers, *_local_to_global_index_map, t, x, M, K, b); } +void HeatConductionProcess::assembleWithJacobianConcreteProcess( + const double t, GlobalVector const& x, GlobalVector const& xdot, + const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K, + GlobalVector& b, GlobalMatrix& Jac) +{ + DBUG("AssembleWithJacobian HeatConductionProcess."); + + // Call global assembler for each local assembly item. + GlobalExecutor::executeMemberDereferenced( + _global_assembler, &VectorMatrixAssembler::assembleWithJacobian, + _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx, + dx_dx, M, K, b, Jac); +} + } // namespace HeatConduction } // namespace ProcessLib diff --git a/ProcessLib/HeatConduction/HeatConductionProcess.h b/ProcessLib/HeatConduction/HeatConductionProcess.h index 280eb5023344a25f60c2cf8e2ba2f79400df710b..d5862b9536a324c7327a7c04d35c05e21dcd939a 100644 --- a/ProcessLib/HeatConduction/HeatConductionProcess.h +++ b/ProcessLib/HeatConduction/HeatConductionProcess.h @@ -25,6 +25,8 @@ class HeatConductionProcess final : public Process public: HeatConductionProcess( MeshLib::Mesh& mesh, + std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& + jacobian_assembler, std::vector<std::unique_ptr<ParameterBase>> const& parameters, std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables, @@ -48,6 +50,11 @@ private: GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) override; + void assembleWithJacobianConcreteProcess( + const double t, GlobalVector const& x, GlobalVector const& xdot, + const double dxdot_dx, const double dx_dx, GlobalMatrix& M, + GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override; + HeatConductionProcessData _process_data; std::vector<std::unique_ptr<HeatConductionLocalAssemblerInterface>>