diff --git a/ProcessLib/ThermoHydroMechanics/Tests.cmake b/ProcessLib/ThermoHydroMechanics/Tests.cmake index 44dc1602d3e5f63ade1b4383b681e50a1a150d4d..a92655671d5a6a36c3cef9a10df4edda53cf9842 100644 --- a/ProcessLib/ThermoHydroMechanics/Tests.cmake +++ b/ProcessLib/ThermoHydroMechanics/Tests.cmake @@ -39,6 +39,24 @@ AddTest( expected_square_1e2_ts_10_t_100.000000.vtu square_1e2_ts_10_t_100.000000.vtu sigma sigma 1e-8 1e-8 ) +# Same as above, but with function instead of group based parameter for Young's modulus +AddTest( +NAME ThermoHydroMechanics_square_1e2_sealed_bimaterial_function +PATH ThermoHydroMechanics/Linear/Beam_sealed_bimaterial +RUNTIME 5 +EXECUTABLE ogs +EXECUTABLE_ARGS square_1e2_function.xml +WRAPPER time +TESTER vtkdiff +REQUIREMENTS NOT OGS_USE_MPI +DIFF_DATA +expected_square_1e2_ts_10_t_100.000000.vtu square_1e2_function_ts_10_t_100.000000.vtu displacement displacement 1e-8 1e-8 +expected_square_1e2_ts_10_t_100.000000.vtu square_1e2_function_ts_10_t_100.000000.vtu pressure pressure 1e-8 1e-8 +expected_square_1e2_ts_10_t_100.000000.vtu square_1e2_function_ts_10_t_100.000000.vtu temperature temperature 1e-8 1e-8 +expected_square_1e2_ts_10_t_100.000000.vtu square_1e2_function_ts_10_t_100.000000.vtu epsilon epsilon 1e-8 1e-8 +expected_square_1e2_ts_10_t_100.000000.vtu square_1e2_function_ts_10_t_100.000000.vtu sigma sigma 1e-8 1e-8 +) + # ThermoHydroMechanics; Small deformation, linear poroelastic, unsealed, bimaterial AddTest( NAME ThermoHydroMechanics_square_1e2_unsealed_bimaterial diff --git a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h index 2206073dfd50d206e82f728665ff21efe07b726f..4a5294b938bdf5b3a178149238042f4752919517 100644 --- a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h +++ b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h @@ -230,9 +230,6 @@ void ThermoHydroMechanicsLocalAssembler< _process_data.solid_materials, _process_data.material_ids, _element.getID()); - ParameterLib::SpatialPosition x_position; - x_position.setElementID(_element.getID()); - auto const& medium = _process_data.media_map->getMedium(_element.getID()); auto const& liquid_phase = medium->phase("AqueousLiquid"); auto const& solid_phase = medium->phase("Solid"); @@ -255,7 +252,6 @@ void ThermoHydroMechanicsLocalAssembler< _integration_method.getNumberOfPoints(); for (unsigned ip = 0; ip < n_integration_points; ip++) { - x_position.setIntegrationPoint(ip); auto const& w = _ip_data[ip].integration_weight; auto const& N_u_op = _ip_data[ip].N_u_op; @@ -273,6 +269,12 @@ void ThermoHydroMechanicsLocalAssembler< auto const T_int_pt = N_T.dot(T); double const dT_int_pt = N_T.dot(T_dot) * dt; + ParameterLib::SpatialPosition const x_position{ + std::nullopt, _element.getID(), ip, + MathLib::Point3d( + NumLib::interpolateCoordinates<ShapeFunctionDisplacement, + ShapeMatricesTypeDisplacement>( + _element, N_u))}; auto const x_coord = NumLib::interpolateXCoordinate<ShapeFunctionDisplacement, ShapeMatricesTypeDisplacement>( @@ -648,9 +650,6 @@ std::vector<double> const& ThermoHydroMechanicsLocalAssembler< temperature_size> const>(local_x.data() + temperature_index, temperature_size); - ParameterLib::SpatialPosition x_position; - x_position.setElementID(_element.getID()); - auto const& medium = _process_data.media_map->getMedium(_element.getID()); auto const& liquid_phase = medium->phase("AqueousLiquid"); auto const& solid_phase = medium->phase("Solid"); @@ -660,10 +659,14 @@ std::vector<double> const& ThermoHydroMechanicsLocalAssembler< for (unsigned ip = 0; ip < n_integration_points; ip++) { - x_position.setIntegrationPoint(ip); - auto const& N_p = _ip_data[ip].N_p; + ParameterLib::SpatialPosition const x_position{ + std::nullopt, _element.getID(), ip, + MathLib::Point3d( + NumLib::interpolateCoordinates<ShapeFunctionDisplacement, + ShapeMatricesTypeDisplacement>( + _element, _ip_data[ip].N_u))}; vars[static_cast<int>(MaterialPropertyLib::Variable::temperature)] = N_p.dot(T); // N_p = N_T double const p_int_pt = N_p.dot(p); @@ -761,8 +764,6 @@ void ThermoHydroMechanicsLocalAssembler< temperature_size> const>(local_xdot.data() + temperature_index, temperature_size); - ParameterLib::SpatialPosition x_position; - x_position.setElementID(_element.getID()); auto const& medium = _process_data.media_map->getMedium(_element.getID()); auto const& solid_phase = medium->phase("Solid"); MaterialPropertyLib::VariableArray vars; @@ -770,11 +771,16 @@ void ThermoHydroMechanicsLocalAssembler< int const n_integration_points = _integration_method.getNumberOfPoints(); for (int ip = 0; ip < n_integration_points; ip++) { - x_position.setIntegrationPoint(ip); auto const& N_u = _ip_data[ip].N_u; auto const& N_T = _ip_data[ip].N_p; auto const& dNdx_u = _ip_data[ip].dNdx_u; + ParameterLib::SpatialPosition const x_position{ + std::nullopt, _element.getID(), ip, + MathLib::Point3d( + NumLib::interpolateCoordinates<ShapeFunctionDisplacement, + ShapeMatricesTypeDisplacement>( + _element, N_u))}; auto const x_coord = NumLib::interpolateXCoordinate<ShapeFunctionDisplacement, ShapeMatricesTypeDisplacement>( diff --git a/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowFEM-impl.h b/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowFEM-impl.h index 66377dbd4c1ec07cb8c3c661a002540974e290b1..73f934e13560120f8b63ed4f177fd22f706b0866 100644 --- a/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowFEM-impl.h +++ b/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowFEM-impl.h @@ -54,12 +54,9 @@ ThermoRichardsFlowLocalAssembler<ShapeFunction, GlobalDim>:: auto const& medium = *_process_data.media_map->getMedium(_element.getID()); - ParameterLib::SpatialPosition x_position; - x_position.setElementID(_element.getID()); for (unsigned ip = 0; ip < n_integration_points; ip++) { auto const& sm = shape_matrices[ip]; - x_position.setIntegrationPoint(ip); _ip_data.emplace_back(); auto& ip_data = _ip_data[ip]; _ip_data[ip].integration_weight = @@ -69,6 +66,11 @@ ThermoRichardsFlowLocalAssembler<ShapeFunction, GlobalDim>:: ip_data.N = sm.N; ip_data.dNdx = sm.dNdx; + ParameterLib::SpatialPosition const x_position{ + std::nullopt, _element.getID(), ip, + MathLib::Point3d(NumLib::interpolateCoordinates<ShapeFunction, + ShapeMatricesType>( + _element, sm.N))}; // Initial porosity. Could be read from integration point data or mesh. ip_data.porosity = medium[MPL::porosity].template initialValue<double>( x_position, @@ -121,17 +123,18 @@ void ThermoRichardsFlowLocalAssembler<ShapeFunction, GlobalDim>:: auto const& medium = *_process_data.media_map->getMedium(_element.getID()); MPL::VariableArray variables; - ParameterLib::SpatialPosition x_position; - x_position.setElementID(_element.getID()); - unsigned const n_integration_points = _integration_method.getNumberOfPoints(); for (unsigned ip = 0; ip < n_integration_points; ip++) { - x_position.setIntegrationPoint(ip); - auto const& N = _ip_data[ip].N; + ParameterLib::SpatialPosition const x_position{ + std::nullopt, _element.getID(), ip, + MathLib::Point3d(NumLib::interpolateCoordinates<ShapeFunction, + ShapeMatricesType>( + _element, N))}; + double p_cap_ip; NumLib::shapeFunctionInterpolate(-p_L, N, p_cap_ip); @@ -225,19 +228,21 @@ void ThermoRichardsFlowLocalAssembler<ShapeFunction, GlobalDim>:: MPL::VariableArray variables; MPL::VariableArray variables_prev; - ParameterLib::SpatialPosition x_position; - x_position.setElementID(_element.getID()); - unsigned const n_integration_points = _integration_method.getNumberOfPoints(); for (unsigned ip = 0; ip < n_integration_points; ip++) { - x_position.setIntegrationPoint(ip); auto const& w = _ip_data[ip].integration_weight; auto const& N = _ip_data[ip].N; auto const& dNdx = _ip_data[ip].dNdx; + ParameterLib::SpatialPosition const x_position{ + std::nullopt, _element.getID(), ip, + MathLib::Point3d(NumLib::interpolateCoordinates<ShapeFunction, + ShapeMatricesType>( + _element, N))}; + double T_ip; NumLib::shapeFunctionInterpolate(T, N, T_ip); double T_dot_ip; @@ -733,19 +738,21 @@ void ThermoRichardsFlowLocalAssembler<ShapeFunction, GlobalDim>::assemble( MPL::VariableArray variables; MPL::VariableArray variables_prev; - ParameterLib::SpatialPosition x_position; - x_position.setElementID(_element.getID()); - unsigned const n_integration_points = _integration_method.getNumberOfPoints(); for (unsigned ip = 0; ip < n_integration_points; ip++) { - x_position.setIntegrationPoint(ip); auto const& w = _ip_data[ip].integration_weight; auto const& N = _ip_data[ip].N; auto const& dNdx = _ip_data[ip].dNdx; + ParameterLib::SpatialPosition const x_position{ + std::nullopt, _element.getID(), ip, + MathLib::Point3d(NumLib::interpolateCoordinates<ShapeFunction, + ShapeMatricesType>( + _element, N))}; + double T_ip; NumLib::shapeFunctionInterpolate(T, N, T_ip); double T_dot_ip; @@ -1231,9 +1238,6 @@ void ThermoRichardsFlowLocalAssembler<ShapeFunction, GlobalDim>:: MPL::VariableArray variables; MPL::VariableArray variables_prev; - ParameterLib::SpatialPosition x_position; - x_position.setElementID(_element.getID()); - unsigned const n_integration_points = _integration_method.getNumberOfPoints(); @@ -1242,9 +1246,14 @@ void ThermoRichardsFlowLocalAssembler<ShapeFunction, GlobalDim>:: for (unsigned ip = 0; ip < n_integration_points; ip++) { - x_position.setIntegrationPoint(ip); auto const& N = _ip_data[ip].N; + ParameterLib::SpatialPosition const x_position{ + std::nullopt, _element.getID(), ip, + MathLib::Point3d(NumLib::interpolateCoordinates<ShapeFunction, + ShapeMatricesType>( + _element, N))}; + double T_ip; NumLib::shapeFunctionInterpolate(T, N, T_ip); double T_dot_ip; diff --git a/ProcessLib/ThermoRichardsMechanics/LocalAssemblerInterface.h b/ProcessLib/ThermoRichardsMechanics/LocalAssemblerInterface.h index 426265183b3371eb7bd609389429c6247a520f71..a13bad27af6452401a6292034902cc5a7bb5dc97 100644 --- a/ProcessLib/ThermoRichardsMechanics/LocalAssemblerInterface.h +++ b/ProcessLib/ThermoRichardsMechanics/LocalAssemblerInterface.h @@ -51,44 +51,6 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface, { material_states_.emplace_back(solid_material_); } - - ParameterLib::SpatialPosition x_position; - x_position.setElementID(e.getID()); - - auto const& medium = process_data_.media_map->getMedium(e.getID()); - - for (unsigned ip = 0; ip < n_integration_points; ip++) - { - namespace MPL = MaterialPropertyLib; - - x_position.setIntegrationPoint(ip); - - auto& current_state = current_states_[ip]; - - // Initial porosity. Could be read from integration point data or - // mesh. - current_state.poro_data.phi = - medium->property(MPL::porosity) - .template initialValue<double>( - x_position, - std::numeric_limits< - double>::quiet_NaN() /* t independent */); - - if (medium->hasProperty(MPL::PropertyType::transport_porosity)) - { - current_state.transport_poro_data.phi = - medium->property(MPL::transport_porosity) - .template initialValue<double>( - x_position, - std::numeric_limits< - double>::quiet_NaN() /* t independent */); - } - else - { - current_state.transport_poro_data.phi = - current_state.poro_data.phi; - } - } } std::size_t setIPDataInitialConditions(std::string const& name, diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h index fa06f3b21868dac80f6f39d5c381a9c0e83e5f23..329301297ce713504be95a96f132e4ff23e74a24 100644 --- a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h +++ b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h @@ -61,11 +61,8 @@ ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunction, DisplacementDim>(e, is_axially_symmetric, integration_method); - ParameterLib::SpatialPosition x_position; - x_position.setElementID(e.getID()); for (unsigned ip = 0; ip < n_integration_points; ip++) { - x_position.setIntegrationPoint(ip); auto& ip_data = ip_data_[ip]; auto const& sm_u = shape_matrices_u[ip]; ip_data_[ip].integration_weight = @@ -117,20 +114,21 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement, this->process_data_.media_map->getMedium(this->element_.getID()); MPL::VariableArray variables; - ParameterLib::SpatialPosition x_position; - x_position.setElementID(this->element_.getID()); - auto const& solid_phase = medium->phase("Solid"); unsigned const n_integration_points = this->integration_method_.getNumberOfPoints(); for (unsigned ip = 0; ip < n_integration_points; ip++) { - x_position.setIntegrationPoint(ip); - // N is used for both T and p variables. auto const& N = ip_data_[ip].N_p; + ParameterLib::SpatialPosition const x_position{ + std::nullopt, this->element_.getID(), ip, + MathLib::Point3d( + NumLib::interpolateCoordinates<ShapeFunctionDisplacement, + ShapeMatricesTypeDisplacement>( + this->element_, ip_data_[ip].N_u))}; double p_cap_ip; NumLib::shapeFunctionInterpolate(-p_L, N, p_cap_ip); @@ -177,9 +175,6 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement, auto& medium = *this->process_data_.media_map->getMedium(this->element_.getID()); - ParameterLib::SpatialPosition x_position; - x_position.setElementID(this->element_.getID()); - LocalMatrices loc_mat; loc_mat.setZero(); LocalMatrices loc_mat_current_ip; @@ -191,7 +186,12 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement, for (unsigned ip = 0; ip < this->integration_method_.getNumberOfPoints(); ++ip) { - x_position.setIntegrationPoint(ip); + ParameterLib::SpatialPosition const x_position{ + std::nullopt, this->element_.getID(), ip, + MathLib::Point3d( + NumLib::interpolateCoordinates<ShapeFunctionDisplacement, + ShapeMatricesTypeDisplacement>( + this->element_, ip_data_[ip].N_u))}; assembleWithJacobianSingleIP( t, dt, x_position, // @@ -463,9 +463,6 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement, auto const& process_data = this->process_data_; auto& medium = *process_data.media_map->getMedium(e_id); - ParameterLib::SpatialPosition x_position; - x_position.setElementID(e_id); - unsigned const n_integration_points = this->integration_method_.getNumberOfPoints(); @@ -490,8 +487,6 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement, auto& current_state = this->current_states_[ip]; auto& output_data = this->output_data_[ip]; - x_position.setIntegrationPoint(ip); - auto const& ip_data = ip_data_[ip]; // N is used for both p and T variables @@ -500,6 +495,12 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement, auto const& dNdx_u = ip_data.dNdx_u; auto const& dNdx = ip_data.dNdx_p; + ParameterLib::SpatialPosition const x_position{ + std::nullopt, this->element_.getID(), ip, + MathLib::Point3d( + NumLib::interpolateCoordinates<ShapeFunctionDisplacement, + ShapeMatricesTypeDisplacement>( + this->element_, N_u))}; auto const x_coord = NumLib::interpolateXCoordinate<ShapeFunctionDisplacement, ShapeMatricesTypeDisplacement>( diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM.h b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM.h index e90d17447464798b97df90ccb7953216f124155b..6166225380c0a4df9dfe621b7c0e0b1d36909ca8 100644 --- a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM.h +++ b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM.h @@ -255,27 +255,51 @@ public: { unsigned const n_integration_points = this->integration_method_.getNumberOfPoints(); + auto const time_independent = std::numeric_limits<double>::quiet_NaN(); + auto const& medium = + *this->process_data_.media_map->getMedium(this->element_.getID()); for (unsigned ip = 0; ip < n_integration_points; ip++) { + ParameterLib::SpatialPosition const x_position{ + std::nullopt, this->element_.getID(), ip, + MathLib::Point3d(NumLib::interpolateCoordinates< + ShapeFunctionDisplacement, + ShapeMatricesTypeDisplacement>( + this->element_, this->ip_data_[ip].N_u))}; + auto& current_state = this->current_states_[ip]; + // Set initial stress from parameter. if (this->process_data_.initial_stress != nullptr) { - ParameterLib::SpatialPosition const x_position{ - std::nullopt, this->element_.getID(), ip, - MathLib::Point3d(NumLib::interpolateCoordinates< - ShapeFunctionDisplacement, - ShapeMatricesTypeDisplacement>( - this->element_, ip_data_[ip].N_u))}; - - this->current_states_[ip].s_mech_data.sigma_eff = + current_state.s_mech_data.sigma_eff = MathLib::KelvinVector::symmetricTensorToKelvinVector< DisplacementDim>((*this->process_data_.initial_stress)( - std::numeric_limits< - double>::quiet_NaN() /* time independent */, - x_position)); + time_independent, x_position)); + } + + // Initial porosity. Could be read from integration point data + // or mesh. + current_state.poro_data.phi = medium.property(MPL::porosity) + .template initialValue<double>( + x_position, time_independent); + + if (medium.hasProperty(MPL::PropertyType::transport_porosity)) + { + current_state.transport_poro_data.phi = + medium.property(MPL::transport_porosity) + .template initialValue<double>(x_position, + time_independent); } + else + { + current_state.transport_poro_data.phi = + current_state.poro_data.phi; + } + } + for (unsigned ip = 0; ip < n_integration_points; ip++) + { this->material_states_[ip].pushBackState(); } diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Beam_sealed_bimaterial/square_1e2_function.xml b/Tests/Data/ThermoHydroMechanics/Linear/Beam_sealed_bimaterial/square_1e2_function.xml new file mode 100644 index 0000000000000000000000000000000000000000..be328afd51d1b02ee3d7417c8fc1e17069de7fc5 --- /dev/null +++ b/Tests/Data/ThermoHydroMechanics/Linear/Beam_sealed_bimaterial/square_1e2_function.xml @@ -0,0 +1,12 @@ +<?xml version='1.0' encoding='ISO-8859-1'?> +<OpenGeoSysProjectDiff base_file="square_1e2.prj"> + <remove sel="/*/parameters/parameter/name[text()="E"]/.." /> + <add sel="/*/parameters"> + <parameter> + <name>E</name> + <type>Function</type> + <expression>if (x < 2) 19; else if (x < 2.1) 0.95; else 9.5;</expression> + </parameter> + </add> + <replace sel="/*/time_loop/output/prefix/text()">square_1e2_function</replace> +</OpenGeoSysProjectDiff> diff --git a/Tests/Data/ThermoHydroMechanics/Linear/ThermoOsmosis/Column.prj b/Tests/Data/ThermoHydroMechanics/Linear/ThermoOsmosis/Column.prj index 68798a516e6223c40e903f5ff14b0b7b8244154e..0e81626b258ad3dac65b432962a89ec62cf88b24 100644 --- a/Tests/Data/ThermoHydroMechanics/Linear/ThermoOsmosis/Column.prj +++ b/Tests/Data/ThermoHydroMechanics/Linear/ThermoOsmosis/Column.prj @@ -217,8 +217,8 @@ <parameters> <parameter> <name>E</name> - <type>Constant</type> - <value>2.88e6</value> + <type>Function</type> + <expression>2.88e6</expression> </parameter> <parameter> <name>nu</name> diff --git a/Tests/Data/ThermoRichardsFlow/ThermoOsmosis/Column.prj b/Tests/Data/ThermoRichardsFlow/ThermoOsmosis/Column.prj index 2d1f0564c727f4cf979ee730953afd10bddf1fae..a90be054612854c159ffd63e0c6285f387eda504 100644 --- a/Tests/Data/ThermoRichardsFlow/ThermoOsmosis/Column.prj +++ b/Tests/Data/ThermoRichardsFlow/ThermoOsmosis/Column.prj @@ -228,8 +228,8 @@ <parameters> <parameter> <name>E</name> - <type>Constant</type> - <value>2.88e6</value> + <type>Function</type> + <expression>2.88e6</expression> </parameter> <parameter> <name>nu</name> diff --git a/Tests/Data/ThermoRichardsMechanics/ThermoOsmosis/Column.prj b/Tests/Data/ThermoRichardsMechanics/ThermoOsmosis/Column.prj index fda993759f275a23baf40de330016d994d6a8379..d20cee32a173e96b2e5e7bba54403fe43ee55b2b 100644 --- a/Tests/Data/ThermoRichardsMechanics/ThermoOsmosis/Column.prj +++ b/Tests/Data/ThermoRichardsMechanics/ThermoOsmosis/Column.prj @@ -233,8 +233,8 @@ <parameters> <parameter> <name>E</name> - <type>Constant</type> - <value>2.88e6</value> + <type>Function</type> + <expression>2.88e6</expression> </parameter> <parameter> <name>nu</name>