From e5520fd77ff0f1df4e9c9a3f7e0db4142b8eb5a0 Mon Sep 17 00:00:00 2001 From: Thomas Fischer Date: Tue, 2 Apr 2024 20:35:44 +0200 Subject: [PATCH 01/21] [PL/LF] Rm templ. parameter in IntegrationPointData Since the shape function isn't stored anymore in IntegrationPointData, the template parameter NodalRowVectorType isn't necessary anymore. --- .../LiquidFlowLocalAssembler-impl.h | 12 ++++------- .../LiquidFlow/LiquidFlowLocalAssembler.h | 21 +++++++------------ 2 files changed, 12 insertions(+), 21 deletions(-) diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h index 4b3e41fffcc..299d76b66e4 100644 --- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h +++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h @@ -333,8 +333,7 @@ void LiquidFlowLocalAssembler::IsotropicCalculator:: calculateLaplacianAndGravityTerm( Eigen::Map& local_K, Eigen::Map& local_b, - IntegrationPointData const& ip_data, + IntegrationPointData const& ip_data, GlobalDimMatrixType const& permeability, double const mu, double const rho_L, GlobalDimVectorType const& specific_body_force, bool const has_gravity) @@ -355,8 +354,7 @@ Eigen::Matrix LiquidFlowLocalAssembler::IsotropicCalculator:: calculateVelocity( Eigen::Map const& local_p, - IntegrationPointData const& ip_data, + IntegrationPointData const& ip_data, GlobalDimMatrixType const& permeability, double const mu, double const rho_L, GlobalDimVectorType const& specific_body_force, bool const has_gravity) @@ -377,8 +375,7 @@ void LiquidFlowLocalAssembler::AnisotropicCalculator:: calculateLaplacianAndGravityTerm( Eigen::Map& local_K, Eigen::Map& local_b, - IntegrationPointData const& ip_data, + IntegrationPointData const& ip_data, GlobalDimMatrixType const& permeability, double const mu, double const rho_L, GlobalDimVectorType const& specific_body_force, bool const has_gravity) @@ -399,8 +396,7 @@ Eigen::Matrix LiquidFlowLocalAssembler::AnisotropicCalculator:: calculateVelocity( Eigen::Map const& local_p, - IntegrationPointData const& ip_data, + IntegrationPointData const& ip_data, GlobalDimMatrixType const& permeability, double const mu, double const rho_L, GlobalDimVectorType const& specific_body_force, bool const has_gravity) diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h index 86ae1838af3..18c361a074a 100644 --- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h +++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h @@ -34,7 +34,7 @@ namespace ProcessLib { namespace LiquidFlow { -template +template struct IntegrationPointData final { explicit IntegrationPointData(GlobalDimNodalMatrixType const& dNdx_, @@ -155,10 +155,9 @@ class LiquidFlowLocalAssembler : public LiquidFlowLocalAssemblerInterface MeshLib::Element const& _element; NumLib::GenericIntegrationMethod const& _integration_method; - std::vector< - IntegrationPointData, - Eigen::aligned_allocator< - IntegrationPointData>> + std::vector, + Eigen::aligned_allocator< + IntegrationPointData>> _ip_data; /** @@ -170,16 +169,14 @@ class LiquidFlowLocalAssembler : public LiquidFlowLocalAssemblerInterface static void calculateLaplacianAndGravityTerm( Eigen::Map& local_K, Eigen::Map& local_b, - IntegrationPointData const& ip_data, + IntegrationPointData const& ip_data, GlobalDimMatrixType const& permeability, double const mu, double const rho_L, GlobalDimVectorType const& specific_body_force, bool const has_gravity); static Eigen::Matrix calculateVelocity( Eigen::Map const& local_p, - IntegrationPointData const& ip_data, + IntegrationPointData const& ip_data, GlobalDimMatrixType const& permeability, double const mu, double const rho_L, GlobalDimVectorType const& specific_body_force, bool const has_gravity); @@ -194,16 +191,14 @@ class LiquidFlowLocalAssembler : public LiquidFlowLocalAssemblerInterface static void calculateLaplacianAndGravityTerm( Eigen::Map& local_K, Eigen::Map& local_b, - IntegrationPointData const& ip_data, + IntegrationPointData const& ip_data, GlobalDimMatrixType const& permeability, double const mu, double const rho_L, GlobalDimVectorType const& specific_body_force, bool const has_gravity); static Eigen::Matrix calculateVelocity( Eigen::Map const& local_p, - IntegrationPointData const& ip_data, + IntegrationPointData const& ip_data, GlobalDimMatrixType const& permeability, double const mu, double const rho_L, GlobalDimVectorType const& specific_body_force, bool const has_gravity); From dac59c6665d0a2af4d0f74d195220ecece7a3816 Mon Sep 17 00:00:00 2001 From: Thomas Fischer Date: Tue, 2 Apr 2024 20:52:28 +0200 Subject: [PATCH 02/21] [PL] Extract shape matrix cache from LiquidFlowProcess to Process --- ProcessLib/LiquidFlow/LiquidFlowProcess.cpp | 1 - ProcessLib/LiquidFlow/LiquidFlowProcess.h | 3 --- ProcessLib/Process.cpp | 1 + ProcessLib/Process.h | 4 ++++ 4 files changed, 5 insertions(+), 4 deletions(-) diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp index 3a5d2af3e8b..1c6fdcfc9d9 100644 --- a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp +++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp @@ -41,7 +41,6 @@ LiquidFlowProcess::LiquidFlowProcess( integration_order, std::move(process_variables), std::move(secondary_variables)), _process_data(std::move(process_data)), - _shape_matrix_cache{integration_order}, _surfaceflux(std::move(surfaceflux)), _is_linear(is_linear) { diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.h b/ProcessLib/LiquidFlow/LiquidFlowProcess.h index aac257a1e6f..8ac5910005e 100644 --- a/ProcessLib/LiquidFlow/LiquidFlowProcess.h +++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.h @@ -17,7 +17,6 @@ #include "LiquidFlowData.h" #include "LiquidFlowLocalAssembler.h" #include "NumLib/DOF/LocalToGlobalIndexMap.h" -#include "NumLib/Fem/ShapeMatrixCache.h" #include "ProcessLib/Process.h" #include "ProcessLib/SurfaceFlux/SurfaceFluxData.h" @@ -109,8 +108,6 @@ class LiquidFlowProcess final : public Process std::vector> _local_assemblers; - NumLib::ShapeMatrixCache _shape_matrix_cache; - std::unique_ptr _surfaceflux; MeshLib::PropertyVector* _hydraulic_flow = nullptr; bool _is_linear = false; diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp index 98ed6566750..1fa91011c3c 100644 --- a/ProcessLib/Process.cpp +++ b/ProcessLib/Process.cpp @@ -71,6 +71,7 @@ Process::Process( } return pcs_BCs; }(_process_variables.size())), + _shape_matrix_cache{integration_order}, _source_term_collections( [&](const std::size_t number_of_processes) -> std::vector diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h index 1816b836668..6086bf62880 100644 --- a/ProcessLib/Process.h +++ b/ProcessLib/Process.h @@ -17,6 +17,7 @@ #include "MaterialLib/MPL/Medium.h" #include "MathLib/LinAlg/GlobalMatrixVectorTypes.h" #include "MeshLib/Utils/IntegrationPointWriter.h" +#include "NumLib/Fem/ShapeMatrixCache.h" #include "NumLib/ODESolver/ODESystem.h" #include "ParameterLib/Parameter.h" #include "ProcessLib/BoundaryConditionAndSourceTerm/BoundaryConditionCollection.h" @@ -394,6 +395,9 @@ class Process /// scheme, the size of vector is the number of the coupled processes. std::vector _boundary_conditions; + /// caches for each mesh element type the shape matrix + NumLib::ShapeMatrixCache _shape_matrix_cache; + private: /// Vector for nodal source term collections. For the monolithic scheme /// or a single process, the size of the vector is one. For the staggered From 405921fc8b7d51a76d688b8ed73a4f03ad09e101 Mon Sep 17 00:00:00 2001 From: Thomas Fischer Date: Tue, 2 Apr 2024 21:12:56 +0200 Subject: [PATCH 03/21] [PL/HT] Pass shape matrix cache to HT local assemblers --- ProcessLib/HT/HTFEM.h | 3 +++ ProcessLib/HT/HTProcess.cpp | 4 ++-- ProcessLib/HT/MonolithicHTFEM.h | 10 ++++++---- ProcessLib/HT/StaggeredHTFEM.h | 9 +++++---- 4 files changed, 16 insertions(+), 10 deletions(-) diff --git a/ProcessLib/HT/HTFEM.h b/ProcessLib/HT/HTFEM.h index da14725acc2..3ed82f8411b 100644 --- a/ProcessLib/HT/HTFEM.h +++ b/ProcessLib/HT/HTFEM.h @@ -50,10 +50,12 @@ class HTFEM : public HTLocalAssemblerInterface NumLib::GenericIntegrationMethod const& integration_method, bool const is_axially_symmetric, HTProcessData const& process_data, + NumLib::ShapeMatrixCache const& shape_matrix_cache, const unsigned dof_per_node) : HTLocalAssemblerInterface(), _element(element), _process_data(process_data), + _shape_matrix_cache(shape_matrix_cache), _integration_method(integration_method) { // This assertion is valid only if all nodal d.o.f. use the same shape @@ -166,6 +168,7 @@ class HTFEM : public HTLocalAssemblerInterface protected: MeshLib::Element const& _element; HTProcessData const& _process_data; + NumLib::ShapeMatrixCache const& _shape_matrix_cache; NumLib::GenericIntegrationMethod const& _integration_method; std::vector< diff --git a/ProcessLib/HT/HTProcess.cpp b/ProcessLib/HT/HTProcess.cpp index f1dd20f1f8b..50971ed954c 100644 --- a/ProcessLib/HT/HTProcess.cpp +++ b/ProcessLib/HT/HTProcess.cpp @@ -56,14 +56,14 @@ void HTProcess::initializeConcreteProcess( ProcessLib::createLocalAssemblers( mesh_space_dimension, mesh.getElements(), dof_table, _local_assemblers, NumLib::IntegrationOrder{integration_order}, - mesh.isAxiallySymmetric(), _process_data); + mesh.isAxiallySymmetric(), _process_data, _shape_matrix_cache); } else { ProcessLib::createLocalAssemblers( mesh_space_dimension, mesh.getElements(), dof_table, _local_assemblers, NumLib::IntegrationOrder{integration_order}, - mesh.isAxiallySymmetric(), _process_data); + mesh.isAxiallySymmetric(), _process_data, _shape_matrix_cache); } _secondary_variables.addSecondaryVariable( diff --git a/ProcessLib/HT/MonolithicHTFEM.h b/ProcessLib/HT/MonolithicHTFEM.h index 902036df8d9..73474c5908e 100644 --- a/ProcessLib/HT/MonolithicHTFEM.h +++ b/ProcessLib/HT/MonolithicHTFEM.h @@ -59,10 +59,12 @@ class MonolithicHTFEM : public HTFEM std::size_t const local_matrix_size, NumLib::GenericIntegrationMethod const& integration_method, bool is_axially_symmetric, - HTProcessData const& process_data) - : HTFEM( - element, local_matrix_size, integration_method, - is_axially_symmetric, process_data, NUM_NODAL_DOF) + HTProcessData const& process_data, + NumLib::ShapeMatrixCache const& shape_matrix_cache) + : HTFEM(element, local_matrix_size, + integration_method, + is_axially_symmetric, process_data, + shape_matrix_cache, NUM_NODAL_DOF) { } diff --git a/ProcessLib/HT/StaggeredHTFEM.h b/ProcessLib/HT/StaggeredHTFEM.h index db88872e2ba..3939dc522b4 100644 --- a/ProcessLib/HT/StaggeredHTFEM.h +++ b/ProcessLib/HT/StaggeredHTFEM.h @@ -58,10 +58,11 @@ class StaggeredHTFEM : public HTFEM std::size_t const local_matrix_size, NumLib::GenericIntegrationMethod const& integration_method, bool is_axially_symmetric, - HTProcessData const& process_data) - : HTFEM(element, local_matrix_size, - integration_method, - is_axially_symmetric, process_data, 1) + HTProcessData const& process_data, + NumLib::ShapeMatrixCache const& shape_matrix_cache) + : HTFEM( + element, local_matrix_size, integration_method, + is_axially_symmetric, process_data, shape_matrix_cache, 1) { } From 2f68e689d3eaad7d2c36aebd19665f2b08990655 Mon Sep 17 00:00:00 2001 From: Thomas Fischer Date: Tue, 2 Apr 2024 21:36:50 +0200 Subject: [PATCH 04/21] [PL/HT] Use shape matrix cache in HT assembler --- ProcessLib/HT/HTFEM.h | 9 +++++++-- ProcessLib/HT/MonolithicHTFEM.h | 6 +++++- ProcessLib/HT/StaggeredHTFEM-impl.h | 12 ++++++++++-- 3 files changed, 22 insertions(+), 5 deletions(-) diff --git a/ProcessLib/HT/HTFEM.h b/ProcessLib/HT/HTFEM.h index 3ed82f8411b..9e2804d3868 100644 --- a/ProcessLib/HT/HTFEM.h +++ b/ProcessLib/HT/HTFEM.h @@ -91,7 +91,8 @@ class HTFEM : public HTLocalAssemblerInterface Eigen::Map getShapeMatrix( const unsigned integration_point) const override { - auto const& N = _ip_data[integration_point].N; + auto const& N = _shape_matrix_cache.NsHigherOrder< + typename ShapeFunction::MeshElement>()[integration_point]; // assumes N is stored contiguously in memory return Eigen::Map(N.data(), N.size()); @@ -273,11 +274,15 @@ class HTFEM : public HTLocalAssemblerInterface *_process_data.media_map.getMedium(_element.getID()); auto const& liquid_phase = medium.phase("AqueousLiquid"); + auto const& Ns = + _shape_matrix_cache + .NsHigherOrder(); + for (unsigned ip = 0; ip < n_integration_points; ++ip) { auto const& ip_data = _ip_data[ip]; - auto const& N = ip_data.N; auto const& dNdx = ip_data.dNdx; + auto const& N = Ns[ip]; pos.setIntegrationPoint(ip); diff --git a/ProcessLib/HT/MonolithicHTFEM.h b/ProcessLib/HT/MonolithicHTFEM.h index 73474c5908e..ff7381ad074 100644 --- a/ProcessLib/HT/MonolithicHTFEM.h +++ b/ProcessLib/HT/MonolithicHTFEM.h @@ -123,13 +123,17 @@ class MonolithicHTFEM : public HTFEM double average_velocity_norm = 0.0; ip_flux_vector.reserve(n_integration_points); + auto const& Ns = + this->_shape_matrix_cache + .template NsHigherOrder(); + for (unsigned ip(0); ip < n_integration_points; ip++) { pos.setIntegrationPoint(ip); auto const& ip_data = this->_ip_data[ip]; - auto const& N = ip_data.N; auto const& dNdx = ip_data.dNdx; + auto const& N = Ns[ip]; auto const& w = ip_data.integration_weight; double T_int_pt = 0.0; diff --git a/ProcessLib/HT/StaggeredHTFEM-impl.h b/ProcessLib/HT/StaggeredHTFEM-impl.h index 65c67a68b1e..7a828cdddac 100644 --- a/ProcessLib/HT/StaggeredHTFEM-impl.h +++ b/ProcessLib/HT/StaggeredHTFEM-impl.h @@ -79,13 +79,17 @@ void StaggeredHTFEM::assembleHydraulicEquation( unsigned const n_integration_points = this->_integration_method.getNumberOfPoints(); + auto const& Ns = + this->_shape_matrix_cache + .template NsHigherOrder(); + for (unsigned ip(0); ip < n_integration_points; ip++) { pos.setIntegrationPoint(ip); auto const& ip_data = this->_ip_data[ip]; - auto const& N = ip_data.N; auto const& dNdx = ip_data.dNdx; + auto const& N = Ns[ip]; auto const& w = ip_data.integration_weight; double p_int_pt = 0.0; @@ -208,13 +212,17 @@ void StaggeredHTFEM::assembleHeatTransportEquation( double average_velocity_norm = 0.0; ip_flux_vector.reserve(n_integration_points); + auto const& Ns = + this->_shape_matrix_cache + .template NsHigherOrder(); + for (unsigned ip(0); ip < n_integration_points; ip++) { pos.setIntegrationPoint(ip); auto const& ip_data = this->_ip_data[ip]; - auto const& N = ip_data.N; auto const& dNdx = ip_data.dNdx; + auto const& N = Ns[ip]; auto const& w = ip_data.integration_weight; double p_at_xi = 0.; From 845cbb1768f564251fd8cd2a1418a18f54f7d390 Mon Sep 17 00:00:00 2001 From: Tom Fischer Date: Wed, 3 Apr 2024 07:27:28 +0200 Subject: [PATCH 05/21] [PL/CT] Pass shape matrix cache to local assemblers of ComponentTransport --- ProcessLib/ComponentTransport/ComponentTransportFEM.h | 4 ++++ ProcessLib/ComponentTransport/ComponentTransportProcess.cpp | 3 ++- 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h index 02ff5f8abff..d60e0c49795 100644 --- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h +++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h @@ -25,6 +25,7 @@ #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h" #include "NumLib/Fem/InitShapeMatrices.h" #include "NumLib/Fem/Integration/GenericIntegrationMethod.h" +#include "NumLib/Fem/ShapeMatrixCache.h" #include "NumLib/Fem/ShapeMatrixPolicy.h" #include "NumLib/Function/Interpolation.h" #include "NumLib/NumericalStability/AdvectionMatrixAssembler.h" @@ -246,6 +247,7 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface NumLib::GenericIntegrationMethod const& integration_method, bool is_axially_symmetric, ComponentTransportProcessData const& process_data, + NumLib::ShapeMatrixCache const& shape_matrix_cache, std::vector> const& transport_process_variables) : temperature_index(process_data.isothermal ? -1 @@ -255,6 +257,7 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface : 2 * ShapeFunction::NPOINTS), _element(element), _process_data(process_data), + _shape_matrix_cache(shape_matrix_cache), _integration_method(integration_method), _transport_process_variables(transport_process_variables) { @@ -1964,6 +1967,7 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface private: MeshLib::Element const& _element; ComponentTransportProcessData const& _process_data; + NumLib::ShapeMatrixCache const& _shape_matrix_cache; NumLib::GenericIntegrationMethod const& _integration_method; std::vector> const diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp index 2e81003a219..dde3556e52b 100644 --- a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp +++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp @@ -23,6 +23,7 @@ #include "MathLib/LinAlg/LinAlg.h" #include "MeshLib/Utils/getOrCreateMeshProperty.h" #include "NumLib/DOF/ComputeSparsityPattern.h" +#include "NumLib/Fem/ShapeMatrixCache.h" #include "NumLib/NumericalStability/FluxCorrectedTransport.h" #include "NumLib/NumericalStability/NumericalStabilization.h" #include "NumLib/ODESolver/NonlinearSystem.h" @@ -165,7 +166,7 @@ void ComponentTransportProcess::initializeConcreteProcess( ProcessLib::createLocalAssemblers( mesh_space_dimension, mesh.getElements(), dof_table, _local_assemblers, NumLib::IntegrationOrder{integration_order}, mesh.isAxiallySymmetric(), - _process_data, transport_process_variables); + _process_data, _shape_matrix_cache, transport_process_variables); if (_chemical_solver_interface && !_use_monolithic_scheme) { From 86b6b483a2dadbfd49b6c7278665a550c6936c39 Mon Sep 17 00:00:00 2001 From: Tom Fischer Date: Wed, 3 Apr 2024 10:11:36 +0200 Subject: [PATCH 06/21] [PL/CT] Use shape matrix cache --- .../ComponentTransportFEM.h | 78 +++++++++++++++---- 1 file changed, 65 insertions(+), 13 deletions(-) diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h index d60e0c49795..c5020a3a9cd 100644 --- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h +++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h @@ -328,12 +328,17 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface ParameterLib::SpatialPosition pos; pos.setElementID(_element.getID()); + auto const& Ns = + _shape_matrix_cache + .NsHigherOrder(); + unsigned const n_integration_points = _integration_method.getNumberOfPoints(); + for (unsigned ip = 0; ip < n_integration_points; ip++) { auto& ip_data = _ip_data[ip]; - auto const& N = ip_data.N; + auto const& N = Ns[ip]; auto const& chemical_system_id = ip_data.chemical_system_id; auto const n_component = _transport_process_variables.size(); @@ -372,12 +377,17 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface ParameterLib::SpatialPosition pos; pos.setElementID(_element.getID()); + auto const& Ns = + _shape_matrix_cache + .NsHigherOrder(); + unsigned const n_integration_points = _integration_method.getNumberOfPoints(); + for (unsigned ip = 0; ip < n_integration_points; ip++) { auto& ip_data = _ip_data[ip]; - auto const& N = ip_data.N; + auto const& N = Ns[ip]; auto& porosity = ip_data.porosity; auto const& porosity_prev = ip_data.porosity_prev; auto const& chemical_system_id = ip_data.chemical_system_id; @@ -595,13 +605,17 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface ip_flux_vector.reserve(n_integration_points); } + auto const& Ns = + _shape_matrix_cache + .NsHigherOrder(); + for (unsigned ip(0); ip < n_integration_points; ++ip) { pos.setIntegrationPoint(ip); auto& ip_data = _ip_data[ip]; - auto const& N = ip_data.N; auto const& dNdx = ip_data.dNdx; + auto const& N = Ns[ip]; auto const& w = ip_data.integration_weight; auto& porosity = ip_data.porosity; @@ -753,11 +767,15 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface auto const& component = phase.component( _transport_process_variables[component_id].get().getName()); + auto const& Ns = + _shape_matrix_cache + .NsHigherOrder(); + for (unsigned ip(0); ip < n_integration_points; ++ip) { auto& ip_data = _ip_data[ip]; - auto const& N = ip_data.N; auto const& w = ip_data.integration_weight; + auto const& N = Ns[ip]; auto& porosity = ip_data.porosity; auto const retardation_factor = @@ -846,14 +864,18 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface MaterialPropertyLib::VariableArray vars; MaterialPropertyLib::VariableArray vars_prev; + auto const& Ns = + _shape_matrix_cache + .NsHigherOrder(); + for (unsigned ip(0); ip < n_integration_points; ++ip) { pos.setIntegrationPoint(ip); auto& ip_data = _ip_data[ip]; - auto const& N = ip_data.N; auto const& dNdx = ip_data.dNdx; auto const& w = ip_data.integration_weight; + auto const& N = Ns[ip]; auto& porosity = ip_data.porosity; auto const& porosity_prev = ip_data.porosity_prev; @@ -970,14 +992,18 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface double average_velocity_norm = 0.0; ip_flux_vector.reserve(n_integration_points); + auto const& Ns = + _shape_matrix_cache + .NsHigherOrder(); + for (unsigned ip(0); ip < n_integration_points; ip++) { pos.setIntegrationPoint(ip); auto const& ip_data = this->_ip_data[ip]; - auto const& N = ip_data.N; auto const& dNdx = ip_data.dNdx; auto const& w = ip_data.integration_weight; + auto const& N = Ns[ip]; double p_at_xi = 0.; NumLib::shapeFunctionInterpolate(local_p, N, p_at_xi); @@ -1112,14 +1138,18 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface auto const& component = phase.component( _transport_process_variables[component_id].get().getName()); + auto const& Ns = + _shape_matrix_cache + .NsHigherOrder(); + for (unsigned ip(0); ip < n_integration_points; ++ip) { pos.setIntegrationPoint(ip); auto& ip_data = _ip_data[ip]; - auto const& N = ip_data.N; auto const& dNdx = ip_data.dNdx; auto const& w = ip_data.integration_weight; + auto const& N = Ns[ip]; auto& porosity = ip_data.porosity; auto const& porosity_prev = ip_data.porosity_prev; @@ -1310,14 +1340,18 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface MaterialPropertyLib::VariableArray vars; MaterialPropertyLib::VariableArray vars_prev; + auto const& Ns = + _shape_matrix_cache + .NsHigherOrder(); + for (unsigned ip(0); ip < n_integration_points; ++ip) { pos.setIntegrationPoint(ip); auto& ip_data = _ip_data[ip]; - auto const& N = ip_data.N; auto const& dNdx = ip_data.dNdx; auto const& w = ip_data.integration_weight; + auto const& N = Ns[ip]; auto& phi = ip_data.porosity; auto const& phi_prev = ip_data.porosity_prev; @@ -1432,14 +1466,18 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface auto const& component = phase.component( _transport_process_variables[component_id].get().getName()); + auto const& Ns = + _shape_matrix_cache + .NsHigherOrder(); + for (unsigned ip(0); ip < n_integration_points; ++ip) { pos.setIntegrationPoint(ip); auto& ip_data = _ip_data[ip]; - auto const& N = ip_data.N; auto const& dNdx = ip_data.dNdx; auto const& w = ip_data.integration_weight; + auto const& N = Ns[ip]; auto& phi = ip_data.porosity; auto const& phi_prev = ip_data.porosity_prev; @@ -1555,13 +1593,18 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface auto const& medium = *_process_data.media_map.getMedium(_element.getID()); auto const component_id = transport_process_id - 1; + + auto const& Ns = + _shape_matrix_cache + .NsHigherOrder(); + for (unsigned ip(0); ip < n_integration_points; ++ip) { pos.setIntegrationPoint(ip); auto& ip_data = _ip_data[ip]; - auto const& N = ip_data.N; auto const w = ip_data.integration_weight; + auto const& N = Ns[ip]; auto& porosity = ip_data.porosity; auto const& porosity_prev = ip_data.porosity_prev; auto const chemical_system_id = ip_data.chemical_system_id; @@ -1672,11 +1715,15 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface *_process_data.media_map.getMedium(_element.getID()); auto const& phase = medium.phase("AqueousLiquid"); + auto const& Ns = + _shape_matrix_cache + .NsHigherOrder(); + for (unsigned ip = 0; ip < n_integration_points; ++ip) { auto const& ip_data = _ip_data[ip]; - auto const& N = ip_data.N; auto const& dNdx = ip_data.dNdx; + auto const& N = Ns[ip]; auto const& porosity = ip_data.porosity; pos.setIntegrationPoint(ip); @@ -1718,7 +1765,8 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface Eigen::Map getShapeMatrix( const unsigned integration_point) const override { - auto const& N = _ip_data[integration_point].N; + auto const& N = _shape_matrix_cache.NsHigherOrder< + typename ShapeFunction::MeshElement>()[integration_point]; // assumes N is stored contiguously in memory return Eigen::Map(N.data(), N.size()); @@ -1899,11 +1947,15 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface auto const& component = phase.component( _transport_process_variables[component_id].get().getName()); + auto const& Ns = + _shape_matrix_cache + .NsHigherOrder(); + for (unsigned ip = 0; ip < n_integration_points; ++ip) { auto const& ip_data = _ip_data[ip]; - auto const& N = ip_data.N; auto const& dNdx = ip_data.dNdx; + auto const& N = Ns[ip]; auto const& phi = ip_data.porosity; pos.setIntegrationPoint(ip); From e88054f381d9bfe415338e3faf3093e3b0fb7673 Mon Sep 17 00:00:00 2001 From: Tom Fischer Date: Wed, 3 Apr 2024 10:35:42 +0200 Subject: [PATCH 07/21] [PL] Pass shape matrix cache to NL::assembleAdvectionMatrix() --- .../AdvectionMatrixAssembler.h | 59 ++++++++++++++++++- .../ComponentTransportFEM.h | 20 ++++--- ProcessLib/HT/MonolithicHTFEM.h | 5 +- ProcessLib/HT/StaggeredHTFEM-impl.h | 5 +- 4 files changed, 75 insertions(+), 14 deletions(-) diff --git a/NumLib/NumericalStability/AdvectionMatrixAssembler.h b/NumLib/NumericalStability/AdvectionMatrixAssembler.h index 3d541b5c048..ed5c09964d7 100644 --- a/NumLib/NumericalStability/AdvectionMatrixAssembler.h +++ b/NumLib/NumericalStability/AdvectionMatrixAssembler.h @@ -14,12 +14,35 @@ #include #include +#include "NumLib/Fem/ShapeMatrixCache.h" #include "NumericalStabilization.h" namespace NumLib { namespace detail { +template +void assembleAdvectionMatrix(IPData const& ip_data_vector, + NumLib::ShapeMatrixCache const& shape_matrix_cache, + std::vector const& ip_flux_vector, + Eigen::MatrixBase& laplacian_matrix) +{ + auto const& Ns = shape_matrix_cache.NsHigherOrder(); + + for (std::size_t ip = 0; ip < ip_flux_vector.size(); ++ip) + { + auto const& ip_data = ip_data_vector[ip]; + auto const w = ip_data.integration_weight; + auto const& dNdx = ip_data.dNdx; + auto const& N = Ns[ip]; + laplacian_matrix.noalias() += + N.transpose() * ip_flux_vector[ip].transpose() * dNdx * w; + } +} + template void assembleAdvectionMatrix(IPData const& ip_data_vector, std::vector const& ip_flux_vector, @@ -79,9 +102,13 @@ void applyFullUpwind(IPData const& ip_data_vector, } } // namespace detail -template +template void assembleAdvectionMatrix(NumericalStabilization const& stabilizer, IPData const& ip_data_vector, + NumLib::ShapeMatrixCache const& shape_matrix_cache, std::vector const& ip_flux_vector, double const average_velocity, Eigen::MatrixBase& laplacian_matrix) @@ -100,8 +127,36 @@ void assembleAdvectionMatrix(NumericalStabilization const& stabilizer, } } - detail::assembleAdvectionMatrix(ip_data_vector, ip_flux_vector, + detail::assembleAdvectionMatrix( + ip_data_vector, shape_matrix_cache, ip_flux_vector, + laplacian_matrix); + }, + stabilizer); +} + +template +void assembleAdvectionMatrix(NumericalStabilization const& stabilizer, + IPData const& ip_data_vector, + std::vector const& ip_flux_vector, + double const average_velocity, + Eigen::MatrixBase& laplacian_matrix) +{ + std::visit( + [&](auto&& stabilizer) + { + using Stabilizer = std::decay_t; + if constexpr (std::is_same_v) + { + if (average_velocity > stabilizer.getCutoffVelocity()) + { + detail::applyFullUpwind(ip_data_vector, ip_flux_vector, laplacian_matrix); + return; + } + } + + detail::assembleAdvectionMatrix( + ip_data_vector, ip_flux_vector, laplacian_matrix); }, stabilizer); } diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h index c5020a3a9cd..4560824aea1 100644 --- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h +++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h @@ -736,9 +736,11 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface if (!_process_data.non_advective_form) { - NumLib::assembleAdvectionMatrix( + NumLib::assembleAdvectionMatrix< + typename ShapeFunction::MeshElement>( _process_data.stabilizer, _ip_data, + _shape_matrix_cache, ip_flux_vector, average_velocity_norm / static_cast(n_integration_points), @@ -1072,8 +1074,9 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface average_velocity_norm += velocity.norm(); } - NumLib::assembleAdvectionMatrix( - process_data.stabilizer, this->_ip_data, ip_flux_vector, + NumLib::assembleAdvectionMatrix( + process_data.stabilizer, this->_ip_data, _shape_matrix_cache, + ip_flux_vector, average_velocity_norm / static_cast(n_integration_points), local_K); } @@ -1273,9 +1276,9 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface if (!_process_data.non_advective_form) { - NumLib::assembleAdvectionMatrix( - _process_data.stabilizer, - _ip_data, + NumLib::assembleAdvectionMatrix< + typename ShapeFunction::MeshElement>( + _process_data.stabilizer, _ip_data, _shape_matrix_cache, ip_flux_vector, average_velocity_norm / static_cast(n_integration_points), @@ -1554,8 +1557,9 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface average_velocity_norm += q.norm(); } - NumLib::assembleAdvectionMatrix( - _process_data.stabilizer, _ip_data, ip_flux_vector, + NumLib::assembleAdvectionMatrix( + _process_data.stabilizer, _ip_data, _shape_matrix_cache, + ip_flux_vector, average_velocity_norm / static_cast(n_integration_points), KCC_Laplacian); diff --git a/ProcessLib/HT/MonolithicHTFEM.h b/ProcessLib/HT/MonolithicHTFEM.h index ff7381ad074..552c69dae28 100644 --- a/ProcessLib/HT/MonolithicHTFEM.h +++ b/ProcessLib/HT/MonolithicHTFEM.h @@ -216,8 +216,9 @@ class MonolithicHTFEM : public HTFEM * in buoyancy effects */ } - NumLib::assembleAdvectionMatrix( - process_data.stabilizer, this->_ip_data, ip_flux_vector, + NumLib::assembleAdvectionMatrix( + process_data.stabilizer, this->_ip_data, this->_shape_matrix_cache, + ip_flux_vector, average_velocity_norm / static_cast(n_integration_points), KTT); } diff --git a/ProcessLib/HT/StaggeredHTFEM-impl.h b/ProcessLib/HT/StaggeredHTFEM-impl.h index 7a828cdddac..d2c5f10c22e 100644 --- a/ProcessLib/HT/StaggeredHTFEM-impl.h +++ b/ProcessLib/HT/StaggeredHTFEM-impl.h @@ -287,8 +287,9 @@ void StaggeredHTFEM::assembleHeatTransportEquation( average_velocity_norm += velocity.norm(); } - NumLib::assembleAdvectionMatrix( - process_data.stabilizer, this->_ip_data, ip_flux_vector, + NumLib::assembleAdvectionMatrix( + process_data.stabilizer, this->_ip_data, this->_shape_matrix_cache, + ip_flux_vector, average_velocity_norm / static_cast(n_integration_points), local_K); } From d949cab8b181ca51c142ef6792a5a9a3b1508208 Mon Sep 17 00:00:00 2001 From: Tom Fischer Date: Wed, 3 Apr 2024 11:09:25 +0200 Subject: [PATCH 08/21] [PL/HT] Rm shape functions per local assembler --- ProcessLib/HT/HTFEM.h | 8 ++------ ProcessLib/HT/HTLocalAssemblerInterface.h | 10 +++------- 2 files changed, 5 insertions(+), 13 deletions(-) diff --git a/ProcessLib/HT/HTFEM.h b/ProcessLib/HT/HTFEM.h index 9e2804d3868..a366e9f6577 100644 --- a/ProcessLib/HT/HTFEM.h +++ b/ProcessLib/HT/HTFEM.h @@ -81,7 +81,7 @@ class HTFEM : public HTLocalAssemblerInterface for (unsigned ip = 0; ip < n_integration_points; ip++) { _ip_data.emplace_back( - shape_matrices[ip].N, shape_matrices[ip].dNdx, + shape_matrices[ip].dNdx, _integration_method.getWeightedPoint(ip).getWeight() * shape_matrices[ip].integralMeasure * shape_matrices[ip].detJ * aperture_size); @@ -172,11 +172,7 @@ class HTFEM : public HTLocalAssemblerInterface NumLib::ShapeMatrixCache const& _shape_matrix_cache; NumLib::GenericIntegrationMethod const& _integration_method; - std::vector< - IntegrationPointData, - Eigen::aligned_allocator< - IntegrationPointData>> - _ip_data; + std::vector> _ip_data; double getHeatEnergyCoefficient( MaterialPropertyLib::VariableArray const& vars, const double porosity, diff --git a/ProcessLib/HT/HTLocalAssemblerInterface.h b/ProcessLib/HT/HTLocalAssemblerInterface.h index c058a23b5c1..8fcd126c7de 100644 --- a/ProcessLib/HT/HTLocalAssemblerInterface.h +++ b/ProcessLib/HT/HTLocalAssemblerInterface.h @@ -21,19 +21,15 @@ namespace ProcessLib namespace HT { -template +template struct IntegrationPointData final { - IntegrationPointData(NodalRowVectorType N_, - GlobalDimNodalMatrixType dNdx_, + IntegrationPointData(GlobalDimNodalMatrixType dNdx_, double const& integration_weight_) - : N(std::move(N_)), - dNdx(std::move(dNdx_)), - integration_weight(integration_weight_) + : dNdx(std::move(dNdx_)), integration_weight(integration_weight_) { } - NodalRowVectorType const N; GlobalDimNodalMatrixType const dNdx; double const integration_weight; From e7c4672f39594d9d230adc9662fb9a193235da2d Mon Sep 17 00:00:00 2001 From: Tom Fischer Date: Wed, 3 Apr 2024 11:22:16 +0200 Subject: [PATCH 09/21] [PL/CT] Rm shape functions per local assembler --- .../ComponentTransport/ComponentTransportFEM.h | 16 +++++----------- 1 file changed, 5 insertions(+), 11 deletions(-) diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h index 4560824aea1..736f81d3aac 100644 --- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h +++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h @@ -38,18 +38,16 @@ namespace ProcessLib { namespace ComponentTransport { -template +template struct IntegrationPointData final { - IntegrationPointData(NodalRowVectorType const& N_, - GlobalDimNodalMatrixType const& dNdx_, + IntegrationPointData(GlobalDimNodalMatrixType const& dNdx_, double const& integration_weight_) - : N(N_), dNdx(dNdx_), integration_weight(integration_weight_) + : dNdx(dNdx_), integration_weight(integration_weight_) { } void pushBackState() { porosity_prev = porosity; } - NodalRowVectorType const N; GlobalDimNodalMatrixType const dNdx; double const integration_weight; @@ -281,7 +279,7 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface for (unsigned ip = 0; ip < n_integration_points; ip++) { _ip_data.emplace_back( - shape_matrices[ip].N, shape_matrices[ip].dNdx, + shape_matrices[ip].dNdx, _integration_method.getWeightedPoint(ip).getWeight() * shape_matrices[ip].integralMeasure * shape_matrices[ip].detJ * aperture_size); @@ -2029,11 +2027,7 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface std::vector> const _transport_process_variables; - std::vector< - IntegrationPointData, - Eigen::aligned_allocator< - IntegrationPointData>> - _ip_data; + std::vector> _ip_data; double getHeatEnergyCoefficient( MaterialPropertyLib::VariableArray const& vars, const double porosity, From 4c1017e6fec953b6959a7155dbf63de26a4e69f2 Mon Sep 17 00:00:00 2001 From: Tom Fischer Date: Fri, 5 Apr 2024 07:50:49 +0200 Subject: [PATCH 10/21] [PL/LF] Quoting allocator isn't necessary anymore --- ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h index 18c361a074a..81593d770da 100644 --- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h +++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h @@ -155,10 +155,7 @@ class LiquidFlowLocalAssembler : public LiquidFlowLocalAssemblerInterface MeshLib::Element const& _element; NumLib::GenericIntegrationMethod const& _integration_method; - std::vector, - Eigen::aligned_allocator< - IntegrationPointData>> - _ip_data; + std::vector> _ip_data; /** * Calculator of the Laplacian and the gravity term for anisotropic From ee7899ed1e52ed9bdf7ce0827487d9e2c98a7dce Mon Sep 17 00:00:00 2001 From: Tom Fischer Date: Fri, 5 Apr 2024 08:05:40 +0200 Subject: [PATCH 11/21] [PL/LF] Add shape matrix cache to process data --- ProcessLib/LiquidFlow/CreateLiquidFlowProcess.cpp | 3 ++- ProcessLib/LiquidFlow/LiquidFlowData.h | 4 ++++ 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/ProcessLib/LiquidFlow/CreateLiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/CreateLiquidFlowProcess.cpp index c25acd8344f..c494c08b93d 100644 --- a/ProcessLib/LiquidFlow/CreateLiquidFlowProcess.cpp +++ b/ProcessLib/LiquidFlow/CreateLiquidFlowProcess.cpp @@ -152,7 +152,8 @@ std::unique_ptr createLiquidFlowProcess( mesh_space_dimension, std::move(specific_body_force), has_gravity, - *aperture_size_parameter}; + *aperture_size_parameter, + NumLib::ShapeMatrixCache{integration_order}}; return std::make_unique( std::move(name), mesh, std::move(jacobian_assembler), parameters, diff --git a/ProcessLib/LiquidFlow/LiquidFlowData.h b/ProcessLib/LiquidFlow/LiquidFlowData.h index b95d1dbad78..b366ccaf2de 100644 --- a/ProcessLib/LiquidFlow/LiquidFlowData.h +++ b/ProcessLib/LiquidFlow/LiquidFlowData.h @@ -13,6 +13,7 @@ #include #include "MaterialLib/MPL/MaterialSpatialDistributionMap.h" +#include "NumLib/Fem/ShapeMatrixCache.h" #include "ParameterLib/Parameter.h" namespace ProcessLib @@ -34,6 +35,9 @@ struct LiquidFlowData final /// It stores aperture size, which is the thickness of 2D element or the /// cross section area of 1D element. For 3D element, the value is set to 1. ParameterLib::Parameter const& aperture_size; + + /// caches for each mesh element type the shape matrix + NumLib::ShapeMatrixCache shape_matrix_cache; }; } // namespace LiquidFlow From deb36bb792049107a0babae3e4cf18ff074e61be Mon Sep 17 00:00:00 2001 From: Tom Fischer Date: Fri, 5 Apr 2024 08:11:02 +0200 Subject: [PATCH 12/21] [PL/LF] Use shape matrix cache from LiquidFlowData --- ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h | 4 ++-- ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h index 299d76b66e4..82404610c26 100644 --- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h +++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h @@ -143,7 +143,7 @@ void LiquidFlowLocalAssembler:: _process_data.element_rotation_matrices[_element.getID()].transpose() * _process_data.specific_body_force; - auto const& N = _shape_matrix_cache + auto const& N = _process_data.shape_matrix_cache .NsHigherOrder(); for (unsigned ip = 0; ip < n_integration_points; ip++) @@ -295,7 +295,7 @@ void LiquidFlowLocalAssembler:: _process_data.element_rotation_matrices[_element.getID()].transpose() * _process_data.specific_body_force; - auto const& N = _shape_matrix_cache + auto const& N = _process_data.shape_matrix_cache .NsHigherOrder(); for (unsigned ip = 0; ip < n_integration_points; ip++) diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h index 81593d770da..096b10eb498 100644 --- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h +++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h @@ -137,7 +137,7 @@ class LiquidFlowLocalAssembler : public LiquidFlowLocalAssemblerInterface const unsigned integration_point) const override { auto const& N = - _shape_matrix_cache + _process_data.shape_matrix_cache .NsHigherOrder(); // assumes N is stored contiguously in memory From d78989d6890db7a764c0a4d839bf9c9f3275035d Mon Sep 17 00:00:00 2001 From: Tom Fischer Date: Fri, 5 Apr 2024 08:15:13 +0200 Subject: [PATCH 13/21] [PL/LF] Use shape matrix cache as in other processes --- .../LiquidFlow/LiquidFlowLocalAssembler-impl.h | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h index 82404610c26..99bc458dfdc 100644 --- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h +++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h @@ -143,15 +143,16 @@ void LiquidFlowLocalAssembler:: _process_data.element_rotation_matrices[_element.getID()].transpose() * _process_data.specific_body_force; - auto const& N = _process_data.shape_matrix_cache - .NsHigherOrder(); + auto const& Ns = _process_data.shape_matrix_cache + .NsHigherOrder(); for (unsigned ip = 0; ip < n_integration_points; ip++) { auto const& ip_data = _ip_data[ip]; + auto const& N = Ns[ip]; double p = 0.; - NumLib::shapeFunctionInterpolate(local_x, N[ip], p); + NumLib::shapeFunctionInterpolate(local_x, N, p); vars.liquid_phase_pressure = p; // Compute density: @@ -176,7 +177,7 @@ void LiquidFlowLocalAssembler:: // Assemble mass matrix, M local_M.noalias() += (porosity * ddensity_dpressure / fluid_density + storage) * - N[ip].transpose() * N[ip] * ip_data.integration_weight; + N.transpose() * N * ip_data.integration_weight; // Compute viscosity: auto const viscosity = @@ -295,14 +296,15 @@ void LiquidFlowLocalAssembler:: _process_data.element_rotation_matrices[_element.getID()].transpose() * _process_data.specific_body_force; - auto const& N = _process_data.shape_matrix_cache - .NsHigherOrder(); + auto const& Ns = _process_data.shape_matrix_cache + .NsHigherOrder(); for (unsigned ip = 0; ip < n_integration_points; ip++) { auto const& ip_data = _ip_data[ip]; + auto const& N = Ns[ip]; double p = 0.; - NumLib::shapeFunctionInterpolate(local_x, N[ip], p); + NumLib::shapeFunctionInterpolate(local_x, N, p); vars.liquid_phase_pressure = p; // Compute density: From 88c0e8c63afa67d81ec3ce68c8d634203213ded9 Mon Sep 17 00:00:00 2001 From: Tom Fischer Date: Fri, 5 Apr 2024 08:30:04 +0200 Subject: [PATCH 14/21] [PL/LF] Rm ref to shape matrix cache from local assemblers --- ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h | 7 ++----- ProcessLib/LiquidFlow/LiquidFlowProcess.cpp | 2 +- 2 files changed, 3 insertions(+), 6 deletions(-) diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h index 096b10eb498..dcdfeb934e4 100644 --- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h +++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h @@ -86,12 +86,10 @@ class LiquidFlowLocalAssembler : public LiquidFlowLocalAssemblerInterface std::size_t const /*local_matrix_size*/, NumLib::GenericIntegrationMethod const& integration_method, bool const is_axially_symmetric, - LiquidFlowData const& process_data, - NumLib::ShapeMatrixCache const& shape_matrix_cache) + LiquidFlowData const& process_data) : _element(element), _integration_method(integration_method), - _process_data(process_data), - _shape_matrix_cache(shape_matrix_cache) + _process_data(process_data) { unsigned const n_integration_points = _integration_method.getNumberOfPoints(); @@ -223,7 +221,6 @@ class LiquidFlowLocalAssembler : public LiquidFlowLocalAssemblerInterface VelocityCacheType& darcy_velocity_at_ips) const; const LiquidFlowData& _process_data; - NumLib::ShapeMatrixCache const& _shape_matrix_cache; }; } // namespace LiquidFlow diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp index 1c6fdcfc9d9..23f7891cc4d 100644 --- a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp +++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp @@ -59,7 +59,7 @@ void LiquidFlowProcess::initializeConcreteProcess( ProcessLib::createLocalAssemblers( mesh_space_dimension, mesh.getElements(), dof_table, _local_assemblers, NumLib::IntegrationOrder{integration_order}, mesh.isAxiallySymmetric(), - _process_data, _shape_matrix_cache); + _process_data); _secondary_variables.addSecondaryVariable( "darcy_velocity", From 104870c72e00de4fe9f987f6aa8250c2652249d6 Mon Sep 17 00:00:00 2001 From: Tom Fischer Date: Fri, 5 Apr 2024 08:21:24 +0200 Subject: [PATCH 15/21] [PL/HT] Add shape matrix cache to process data --- ProcessLib/HT/CreateHTProcess.cpp | 3 ++- ProcessLib/HT/HTProcessData.h | 4 ++++ 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/ProcessLib/HT/CreateHTProcess.cpp b/ProcessLib/HT/CreateHTProcess.cpp index f8cc3241979..8a5a4b6d44d 100644 --- a/ProcessLib/HT/CreateHTProcess.cpp +++ b/ProcessLib/HT/CreateHTProcess.cpp @@ -218,7 +218,8 @@ std::unique_ptr createHTProcess( std::move(stabilizer), projected_specific_body_force_vectors, mesh_space_dimension, - *aperture_size_parameter}; + *aperture_size_parameter, + NumLib::ShapeMatrixCache{integration_order}}; SecondaryVariableCollection secondary_variables; diff --git a/ProcessLib/HT/HTProcessData.h b/ProcessLib/HT/HTProcessData.h index f12e2c81c79..c951cf6b989 100644 --- a/ProcessLib/HT/HTProcessData.h +++ b/ProcessLib/HT/HTProcessData.h @@ -14,6 +14,7 @@ #include #include "MaterialLib/MPL/MaterialSpatialDistributionMap.h" +#include "NumLib/Fem/ShapeMatrixCache.h" #include "NumLib/NumericalStability/NumericalStabilization.h" #include "ParameterLib/ConstantParameter.h" #include "ParameterLib/Parameter.h" @@ -43,6 +44,9 @@ struct HTProcessData final /// cross section area of 1D element. For 3D element, the value is set to 1. ParameterLib::Parameter const& aperture_size = ParameterLib::ConstantParameter("constant_one", 1.0); + + /// caches for each mesh element type the shape matrix + NumLib::ShapeMatrixCache shape_matrix_cache; }; } // namespace HT From 9a6144c5188b2182bff8a5c1bdfd4ed04fc68e38 Mon Sep 17 00:00:00 2001 From: Tom Fischer Date: Fri, 5 Apr 2024 08:51:25 +0200 Subject: [PATCH 16/21] [PL/HT] Use shape matrix cache from HTProcessData --- ProcessLib/HT/HTFEM.h | 2 +- ProcessLib/HT/MonolithicHTFEM.h | 6 +++--- ProcessLib/HT/StaggeredHTFEM-impl.h | 8 ++++---- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/ProcessLib/HT/HTFEM.h b/ProcessLib/HT/HTFEM.h index a366e9f6577..ca9c204198c 100644 --- a/ProcessLib/HT/HTFEM.h +++ b/ProcessLib/HT/HTFEM.h @@ -91,7 +91,7 @@ class HTFEM : public HTLocalAssemblerInterface Eigen::Map getShapeMatrix( const unsigned integration_point) const override { - auto const& N = _shape_matrix_cache.NsHigherOrder< + auto const& N = _process_data.shape_matrix_cache.NsHigherOrder< typename ShapeFunction::MeshElement>()[integration_point]; // assumes N is stored contiguously in memory diff --git a/ProcessLib/HT/MonolithicHTFEM.h b/ProcessLib/HT/MonolithicHTFEM.h index 552c69dae28..d1c741298c6 100644 --- a/ProcessLib/HT/MonolithicHTFEM.h +++ b/ProcessLib/HT/MonolithicHTFEM.h @@ -124,7 +124,7 @@ class MonolithicHTFEM : public HTFEM ip_flux_vector.reserve(n_integration_points); auto const& Ns = - this->_shape_matrix_cache + process_data.shape_matrix_cache .template NsHigherOrder(); for (unsigned ip(0); ip < n_integration_points; ip++) @@ -217,8 +217,8 @@ class MonolithicHTFEM : public HTFEM } NumLib::assembleAdvectionMatrix( - process_data.stabilizer, this->_ip_data, this->_shape_matrix_cache, - ip_flux_vector, + process_data.stabilizer, this->_ip_data, + process_data.shape_matrix_cache, ip_flux_vector, average_velocity_norm / static_cast(n_integration_points), KTT); } diff --git a/ProcessLib/HT/StaggeredHTFEM-impl.h b/ProcessLib/HT/StaggeredHTFEM-impl.h index d2c5f10c22e..7f92bddc13b 100644 --- a/ProcessLib/HT/StaggeredHTFEM-impl.h +++ b/ProcessLib/HT/StaggeredHTFEM-impl.h @@ -80,7 +80,7 @@ void StaggeredHTFEM::assembleHydraulicEquation( this->_integration_method.getNumberOfPoints(); auto const& Ns = - this->_shape_matrix_cache + process_data.shape_matrix_cache .template NsHigherOrder(); for (unsigned ip(0); ip < n_integration_points; ip++) @@ -213,7 +213,7 @@ void StaggeredHTFEM::assembleHeatTransportEquation( ip_flux_vector.reserve(n_integration_points); auto const& Ns = - this->_shape_matrix_cache + process_data.shape_matrix_cache .template NsHigherOrder(); for (unsigned ip(0); ip < n_integration_points; ip++) @@ -288,8 +288,8 @@ void StaggeredHTFEM::assembleHeatTransportEquation( } NumLib::assembleAdvectionMatrix( - process_data.stabilizer, this->_ip_data, this->_shape_matrix_cache, - ip_flux_vector, + process_data.stabilizer, this->_ip_data, + process_data.shape_matrix_cache, ip_flux_vector, average_velocity_norm / static_cast(n_integration_points), local_K); } From b68eee2c56607a9ebfb277963b023d679fe688eb Mon Sep 17 00:00:00 2001 From: Tom Fischer Date: Fri, 5 Apr 2024 08:52:42 +0200 Subject: [PATCH 17/21] [PL/HT] Rm ref to shape matrix cache from local assemblers --- ProcessLib/HT/HTFEM.h | 5 +---- ProcessLib/HT/HTProcess.cpp | 4 ++-- ProcessLib/HT/MonolithicHTFEM.h | 10 ++++------ ProcessLib/HT/StaggeredHTFEM.h | 9 ++++----- 4 files changed, 11 insertions(+), 17 deletions(-) diff --git a/ProcessLib/HT/HTFEM.h b/ProcessLib/HT/HTFEM.h index ca9c204198c..f0e08207fa6 100644 --- a/ProcessLib/HT/HTFEM.h +++ b/ProcessLib/HT/HTFEM.h @@ -50,12 +50,10 @@ class HTFEM : public HTLocalAssemblerInterface NumLib::GenericIntegrationMethod const& integration_method, bool const is_axially_symmetric, HTProcessData const& process_data, - NumLib::ShapeMatrixCache const& shape_matrix_cache, const unsigned dof_per_node) : HTLocalAssemblerInterface(), _element(element), _process_data(process_data), - _shape_matrix_cache(shape_matrix_cache), _integration_method(integration_method) { // This assertion is valid only if all nodal d.o.f. use the same shape @@ -169,7 +167,6 @@ class HTFEM : public HTLocalAssemblerInterface protected: MeshLib::Element const& _element; HTProcessData const& _process_data; - NumLib::ShapeMatrixCache const& _shape_matrix_cache; NumLib::GenericIntegrationMethod const& _integration_method; std::vector> _ip_data; @@ -271,7 +268,7 @@ class HTFEM : public HTLocalAssemblerInterface auto const& liquid_phase = medium.phase("AqueousLiquid"); auto const& Ns = - _shape_matrix_cache + _process_data.shape_matrix_cache .NsHigherOrder(); for (unsigned ip = 0; ip < n_integration_points; ++ip) diff --git a/ProcessLib/HT/HTProcess.cpp b/ProcessLib/HT/HTProcess.cpp index 50971ed954c..f1dd20f1f8b 100644 --- a/ProcessLib/HT/HTProcess.cpp +++ b/ProcessLib/HT/HTProcess.cpp @@ -56,14 +56,14 @@ void HTProcess::initializeConcreteProcess( ProcessLib::createLocalAssemblers( mesh_space_dimension, mesh.getElements(), dof_table, _local_assemblers, NumLib::IntegrationOrder{integration_order}, - mesh.isAxiallySymmetric(), _process_data, _shape_matrix_cache); + mesh.isAxiallySymmetric(), _process_data); } else { ProcessLib::createLocalAssemblers( mesh_space_dimension, mesh.getElements(), dof_table, _local_assemblers, NumLib::IntegrationOrder{integration_order}, - mesh.isAxiallySymmetric(), _process_data, _shape_matrix_cache); + mesh.isAxiallySymmetric(), _process_data); } _secondary_variables.addSecondaryVariable( diff --git a/ProcessLib/HT/MonolithicHTFEM.h b/ProcessLib/HT/MonolithicHTFEM.h index d1c741298c6..91ad8249837 100644 --- a/ProcessLib/HT/MonolithicHTFEM.h +++ b/ProcessLib/HT/MonolithicHTFEM.h @@ -59,12 +59,10 @@ class MonolithicHTFEM : public HTFEM std::size_t const local_matrix_size, NumLib::GenericIntegrationMethod const& integration_method, bool is_axially_symmetric, - HTProcessData const& process_data, - NumLib::ShapeMatrixCache const& shape_matrix_cache) - : HTFEM(element, local_matrix_size, - integration_method, - is_axially_symmetric, process_data, - shape_matrix_cache, NUM_NODAL_DOF) + HTProcessData const& process_data) + : HTFEM( + element, local_matrix_size, integration_method, + is_axially_symmetric, process_data, NUM_NODAL_DOF) { } diff --git a/ProcessLib/HT/StaggeredHTFEM.h b/ProcessLib/HT/StaggeredHTFEM.h index 3939dc522b4..db88872e2ba 100644 --- a/ProcessLib/HT/StaggeredHTFEM.h +++ b/ProcessLib/HT/StaggeredHTFEM.h @@ -58,11 +58,10 @@ class StaggeredHTFEM : public HTFEM std::size_t const local_matrix_size, NumLib::GenericIntegrationMethod const& integration_method, bool is_axially_symmetric, - HTProcessData const& process_data, - NumLib::ShapeMatrixCache const& shape_matrix_cache) - : HTFEM( - element, local_matrix_size, integration_method, - is_axially_symmetric, process_data, shape_matrix_cache, 1) + HTProcessData const& process_data) + : HTFEM(element, local_matrix_size, + integration_method, + is_axially_symmetric, process_data, 1) { } From 3df5356cfafc7c3ad5f84cd2ce376c460c68391a Mon Sep 17 00:00:00 2001 From: Tom Fischer Date: Fri, 5 Apr 2024 09:12:00 +0200 Subject: [PATCH 18/21] [PL/CT] Add shape matrix cache to process data --- ProcessLib/ComponentTransport/ComponentTransportProcessData.h | 4 ++++ .../ComponentTransport/CreateComponentTransportProcess.cpp | 2 +- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcessData.h b/ProcessLib/ComponentTransport/ComponentTransportProcessData.h index 2ac5512eed9..121c3e7641c 100644 --- a/ProcessLib/ComponentTransport/ComponentTransportProcessData.h +++ b/ProcessLib/ComponentTransport/ComponentTransportProcessData.h @@ -17,6 +17,7 @@ #include "LookupTable.h" #include "MaterialLib/MPL/MaterialSpatialDistributionMap.h" #include "MathLib/LinAlg/Eigen/EigenMapTools.h" +#include "NumLib/Fem/ShapeMatrixCache.h" #include "NumLib/NumericalStability/NumericalStabilization.h" #include "ParameterLib/ConstantParameter.h" #include "ParameterLib/Parameter.h" @@ -84,6 +85,9 @@ struct ComponentTransportProcessData bool const isothermal; + /// caches for each mesh element type the shape matrix + NumLib::ShapeMatrixCache shape_matrix_cache; + static const int hydraulic_process_id = 0; // Thermal process is optional, indicated by -1. If present, it is positive. const int thermal_process_id = isothermal ? -1 : 1; diff --git a/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp index bb69b4fe94e..5710f1e4d51 100644 --- a/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp +++ b/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp @@ -317,7 +317,7 @@ std::unique_ptr createComponentTransportProcess( mesh_space_dimension, *aperture_size_parameter, isothermal, - }; + NumLib::ShapeMatrixCache{integration_order}}; SecondaryVariableCollection secondary_variables; From 4c0404ff09799a651c6888bd80e90596337bba4c Mon Sep 17 00:00:00 2001 From: Tom Fischer Date: Fri, 5 Apr 2024 09:13:01 +0200 Subject: [PATCH 19/21] [PL/CT] Use shape matrix cache from ComponentTransportProcessData --- .../ComponentTransportFEM.h | 40 +++++++++---------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h index 736f81d3aac..c0a44f5bc83 100644 --- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h +++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h @@ -327,7 +327,7 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface pos.setElementID(_element.getID()); auto const& Ns = - _shape_matrix_cache + _process_data.shape_matrix_cache .NsHigherOrder(); unsigned const n_integration_points = @@ -376,7 +376,7 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface pos.setElementID(_element.getID()); auto const& Ns = - _shape_matrix_cache + _process_data.shape_matrix_cache .NsHigherOrder(); unsigned const n_integration_points = @@ -604,7 +604,7 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface } auto const& Ns = - _shape_matrix_cache + _process_data.shape_matrix_cache .NsHigherOrder(); for (unsigned ip(0); ip < n_integration_points; ++ip) @@ -738,7 +738,7 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface typename ShapeFunction::MeshElement>( _process_data.stabilizer, _ip_data, - _shape_matrix_cache, + _process_data.shape_matrix_cache, ip_flux_vector, average_velocity_norm / static_cast(n_integration_points), @@ -768,7 +768,7 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface _transport_process_variables[component_id].get().getName()); auto const& Ns = - _shape_matrix_cache + _process_data.shape_matrix_cache .NsHigherOrder(); for (unsigned ip(0); ip < n_integration_points; ++ip) @@ -865,7 +865,7 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface MaterialPropertyLib::VariableArray vars_prev; auto const& Ns = - _shape_matrix_cache + _process_data.shape_matrix_cache .NsHigherOrder(); for (unsigned ip(0); ip < n_integration_points; ++ip) @@ -993,7 +993,7 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface ip_flux_vector.reserve(n_integration_points); auto const& Ns = - _shape_matrix_cache + _process_data.shape_matrix_cache .NsHigherOrder(); for (unsigned ip(0); ip < n_integration_points; ip++) @@ -1073,8 +1073,8 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface } NumLib::assembleAdvectionMatrix( - process_data.stabilizer, this->_ip_data, _shape_matrix_cache, - ip_flux_vector, + process_data.stabilizer, this->_ip_data, + _process_data.shape_matrix_cache, ip_flux_vector, average_velocity_norm / static_cast(n_integration_points), local_K); } @@ -1140,7 +1140,7 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface _transport_process_variables[component_id].get().getName()); auto const& Ns = - _shape_matrix_cache + _process_data.shape_matrix_cache .NsHigherOrder(); for (unsigned ip(0); ip < n_integration_points; ++ip) @@ -1276,8 +1276,8 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface { NumLib::assembleAdvectionMatrix< typename ShapeFunction::MeshElement>( - _process_data.stabilizer, _ip_data, _shape_matrix_cache, - ip_flux_vector, + _process_data.stabilizer, _ip_data, + _process_data.shape_matrix_cache, ip_flux_vector, average_velocity_norm / static_cast(n_integration_points), KCC_Laplacian); @@ -1342,7 +1342,7 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface MaterialPropertyLib::VariableArray vars_prev; auto const& Ns = - _shape_matrix_cache + _process_data.shape_matrix_cache .NsHigherOrder(); for (unsigned ip(0); ip < n_integration_points; ++ip) @@ -1468,7 +1468,7 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface _transport_process_variables[component_id].get().getName()); auto const& Ns = - _shape_matrix_cache + _process_data.shape_matrix_cache .NsHigherOrder(); for (unsigned ip(0); ip < n_integration_points; ++ip) @@ -1556,8 +1556,8 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface } NumLib::assembleAdvectionMatrix( - _process_data.stabilizer, _ip_data, _shape_matrix_cache, - ip_flux_vector, + _process_data.stabilizer, _ip_data, + _process_data.shape_matrix_cache, ip_flux_vector, average_velocity_norm / static_cast(n_integration_points), KCC_Laplacian); @@ -1597,7 +1597,7 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface auto const component_id = transport_process_id - 1; auto const& Ns = - _shape_matrix_cache + _process_data.shape_matrix_cache .NsHigherOrder(); for (unsigned ip(0); ip < n_integration_points; ++ip) @@ -1718,7 +1718,7 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface auto const& phase = medium.phase("AqueousLiquid"); auto const& Ns = - _shape_matrix_cache + _process_data.shape_matrix_cache .NsHigherOrder(); for (unsigned ip = 0; ip < n_integration_points; ++ip) @@ -1767,7 +1767,7 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface Eigen::Map getShapeMatrix( const unsigned integration_point) const override { - auto const& N = _shape_matrix_cache.NsHigherOrder< + auto const& N = _process_data.shape_matrix_cache.NsHigherOrder< typename ShapeFunction::MeshElement>()[integration_point]; // assumes N is stored contiguously in memory @@ -1950,7 +1950,7 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface _transport_process_variables[component_id].get().getName()); auto const& Ns = - _shape_matrix_cache + _process_data.shape_matrix_cache .NsHigherOrder(); for (unsigned ip = 0; ip < n_integration_points; ++ip) From a8550fe0e45340c23d964659af7a718c25c0dfc2 Mon Sep 17 00:00:00 2001 From: Tom Fischer Date: Fri, 5 Apr 2024 09:17:04 +0200 Subject: [PATCH 20/21] [PL/CT] Rm ref to shape matrix cache from local assemblers --- ProcessLib/ComponentTransport/ComponentTransportFEM.h | 3 --- ProcessLib/ComponentTransport/ComponentTransportProcess.cpp | 2 +- 2 files changed, 1 insertion(+), 4 deletions(-) diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h index c0a44f5bc83..d9cedbbc2f5 100644 --- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h +++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h @@ -245,7 +245,6 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface NumLib::GenericIntegrationMethod const& integration_method, bool is_axially_symmetric, ComponentTransportProcessData const& process_data, - NumLib::ShapeMatrixCache const& shape_matrix_cache, std::vector> const& transport_process_variables) : temperature_index(process_data.isothermal ? -1 @@ -255,7 +254,6 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface : 2 * ShapeFunction::NPOINTS), _element(element), _process_data(process_data), - _shape_matrix_cache(shape_matrix_cache), _integration_method(integration_method), _transport_process_variables(transport_process_variables) { @@ -2021,7 +2019,6 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface private: MeshLib::Element const& _element; ComponentTransportProcessData const& _process_data; - NumLib::ShapeMatrixCache const& _shape_matrix_cache; NumLib::GenericIntegrationMethod const& _integration_method; std::vector> const diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp index dde3556e52b..74cfee620f7 100644 --- a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp +++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp @@ -166,7 +166,7 @@ void ComponentTransportProcess::initializeConcreteProcess( ProcessLib::createLocalAssemblers( mesh_space_dimension, mesh.getElements(), dof_table, _local_assemblers, NumLib::IntegrationOrder{integration_order}, mesh.isAxiallySymmetric(), - _process_data, _shape_matrix_cache, transport_process_variables); + _process_data, transport_process_variables); if (_chemical_solver_interface && !_use_monolithic_scheme) { From 64cf5e72fe851ae21992791f7b45c5629b6f10c0 Mon Sep 17 00:00:00 2001 From: Tom Fischer Date: Fri, 5 Apr 2024 09:26:05 +0200 Subject: [PATCH 21/21] [PL] Rm shape matrix cache instance from Process shape matrix cache is now stored in the ProcessData of the concrete process --- ProcessLib/Process.cpp | 1 - ProcessLib/Process.h | 4 ---- 2 files changed, 5 deletions(-) diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp index 1fa91011c3c..98ed6566750 100644 --- a/ProcessLib/Process.cpp +++ b/ProcessLib/Process.cpp @@ -71,7 +71,6 @@ Process::Process( } return pcs_BCs; }(_process_variables.size())), - _shape_matrix_cache{integration_order}, _source_term_collections( [&](const std::size_t number_of_processes) -> std::vector diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h index 6086bf62880..1816b836668 100644 --- a/ProcessLib/Process.h +++ b/ProcessLib/Process.h @@ -17,7 +17,6 @@ #include "MaterialLib/MPL/Medium.h" #include "MathLib/LinAlg/GlobalMatrixVectorTypes.h" #include "MeshLib/Utils/IntegrationPointWriter.h" -#include "NumLib/Fem/ShapeMatrixCache.h" #include "NumLib/ODESolver/ODESystem.h" #include "ParameterLib/Parameter.h" #include "ProcessLib/BoundaryConditionAndSourceTerm/BoundaryConditionCollection.h" @@ -395,9 +394,6 @@ class Process /// scheme, the size of vector is the number of the coupled processes. std::vector _boundary_conditions; - /// caches for each mesh element type the shape matrix - NumLib::ShapeMatrixCache _shape_matrix_cache; - private: /// Vector for nodal source term collections. For the monolithic scheme /// or a single process, the size of the vector is one. For the staggered