Commit 6c00b461 authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#1352 Fix MixedSolidIncompressibility operator, which was not properly updated in #1189.

parent 354558f8
......@@ -183,7 +183,6 @@ namespace MoReFEM
* \param[in] test_pressure_ref_felt Test reference finite element for pressure.
*/
void PrepareInternalDataForQuadraturePoint(const InformationsAtQuadraturePoint& infos_at_quad_pt,
const InformationsAtQuadraturePoint& infos_at_quad_pt_for_test_unknown,
const GeometricElt& geom_elt,
const Advanced::RefFEltInLocalOperator& displacement_ref_felt,
const Advanced::RefFEltInLocalOperator& pressure_ref_felt,
......@@ -218,7 +217,6 @@ namespace MoReFEM
*/
template<unsigned int FeltSpaceDimensionT>
void PrepareInternalDataForQuadraturePointForDimension(const InformationsAtQuadraturePoint& infos_at_quad_pt,
const InformationsAtQuadraturePoint& infos_at_quad_pt_for_test_unknown,
const GeometricElt& geom_elt,
const Advanced::RefFEltInLocalOperator& displacement_ref_felt,
const Advanced::RefFEltInLocalOperator& pressure_ref_felt,
......
......@@ -186,8 +186,6 @@ namespace MoReFEM
auto& tangent_vector_pres_disp = this->vector_parent::template GetLocalVector<EnumUnderlyingType(LocalVectorIndex::tangent_vector_pres_disp)>();
const auto& infos_at_quad_pt_list = elementary_data.GetInformationsAtQuadraturePointList();
assert(infos_at_quad_pt_list.size() == infos_at_quad_pt_list_for_test_unknown.size());
const auto& displacement_ref_felt = elementary_data.GetRefFElt(GetNthUnknown(0));
const auto& pressure_ref_felt = elementary_data.GetRefFElt(GetNthUnknown(1));
......@@ -199,18 +197,13 @@ namespace MoReFEM
const auto felt_space_dimension = displacement_ref_felt.GetFEltSpaceDimension();
unsigned int quad_pt_index = 0;
for (const auto& infos_at_quad_pt : infos_at_quad_pt_list)
{
const auto& infos_at_quad_pt_for_test_unknown = infos_at_quad_pt_list_for_test_unknown[quad_pt_index];
tangent_matrix_disp_disp.Zero();
tangent_vector_disp_pres.Zero();
tangent_vector_pres_disp.Zero();
PrepareInternalDataForQuadraturePoint(infos_at_quad_pt,
infos_at_quad_pt_for_test_unknown,
geom_elt,
displacement_ref_felt,
pressure_ref_felt,
......@@ -219,16 +212,14 @@ namespace MoReFEM
local_displacement,
local_pressure,
felt_space_dimension);
++quad_pt_index;
}
}
template<class HydrostaticLawPolicyT>
void MixedSolidIncompressibility<HydrostaticLawPolicyT>
::PrepareInternalDataForQuadraturePoint(const InformationsAtQuadraturePoint& infos_at_quad_pt,
const InformationsAtQuadraturePoint& infos_at_quad_pt_for_test_unknown,
const GeometricElt& geom_elt,
const Advanced::RefFEltInLocalOperator& displacement_ref_felt,
const Advanced::RefFEltInLocalOperator& pressure_ref_felt,
......@@ -243,7 +234,6 @@ namespace MoReFEM
case 2:
{
this->PrepareInternalDataForQuadraturePointForDimension<2>(infos_at_quad_pt,
infos_at_quad_pt_for_test_unknown,
geom_elt,
displacement_ref_felt,
pressure_ref_felt,
......@@ -256,7 +246,6 @@ namespace MoReFEM
case 3:
{
this->PrepareInternalDataForQuadraturePointForDimension<3>(infos_at_quad_pt,
infos_at_quad_pt_for_test_unknown,
geom_elt,
displacement_ref_felt,
pressure_ref_felt,
......@@ -275,7 +264,6 @@ namespace MoReFEM
template<unsigned int FeltSpaceDimensionT>
void MixedSolidIncompressibility<HydrostaticLawPolicyT>
::PrepareInternalDataForQuadraturePointForDimension(const InformationsAtQuadraturePoint& infos_at_quad_pt,
const InformationsAtQuadraturePoint& infos_at_quad_pt_for_test_unknown,
const GeometricElt& geom_elt,
const Advanced::RefFEltInLocalOperator& displacement_ref_felt,
const Advanced::RefFEltInLocalOperator& pressure_ref_felt,
......@@ -303,8 +291,11 @@ namespace MoReFEM
const int int_Nnode_test_pressure = static_cast<int>(Nnode_test_pressure);
const auto& quad_pt = infos_at_quad_pt.GetQuadraturePoint();
decltype(auto) quad_pt_unknown_list_data = infos_at_quad_pt.GetUnknownData();
decltype(auto) quad_pt_test_unknown_list_data = infos_at_quad_pt.GetTestUnknownData();
const auto weight_meas = quad_pt.GetWeight() * infos_at_quad_pt.GetAbsoluteValueJacobianDeterminant();
const auto weight_meas = quad_pt.GetWeight() * quad_pt_unknown_list_data.GetAbsoluteValueJacobianDeterminant();
auto& tangent_matrix_disp_disp = this->matrix_parent::template GetLocalMatrix<EnumUnderlyingType(LocalMatrixIndex::tangent_matrix_disp_disp)>();
......@@ -344,22 +335,22 @@ namespace MoReFEM
auto& block_contribution_pres_disp = this->matrix_parent::template
GetLocalMatrix<EnumUnderlyingType(LocalMatrixIndex::block_contribution_pres_disp)>();
const auto& grad_felt_phi = infos_at_quad_pt.GetGradientFEltPhi(); // on (u p)
const auto& test_grad_felt_phi = test_unknown_data.GetGradientFEltPhi(); // on (u* p*)
const auto& grad_felt_phi = quad_pt_unknown_list_data.GetGradientFEltPhi(); // on (u p)
const auto& test_grad_felt_phi = quad_pt_test_unknown_list_data.GetGradientFEltPhi(); // on (u* p*)
const auto& dphi_displacement = ExtractSubMatrix(grad_felt_phi, displacement_ref_felt);
const auto& dphi_test_displacement = ExtractSubMatrix(test_grad_felt_phi, test_displacement_ref_felt);
const auto& phi = infos_at_quad_pt.GetFEltPhi(); // on (u p)
const auto& test_phi = infos_at_quad_pt.GetFEltPhi(); // on (u* p*)
const auto& phi = quad_pt_unknown_list_data.GetFEltPhi(); // on (u p)
const auto& test_phi = quad_pt_test_unknown_list_data.GetFEltPhi(); // on (u* p*)
const auto& pressure_phi = ExtractSubVector(phi, pressure_ref_felt);
const auto& test_pressure_phi = ExtractSubVector(test_phi, test_pressure_ref_felt);
Wrappers::Seldon::Transpose(dphi_displacement, transposed_dphi_displacement);
Advanced::OperatorNS::ComputeGradientDisplacementMatrix(infos_at_quad_pt,
Advanced::OperatorNS::ComputeGradientDisplacementMatrix(quad_pt_unknown_list_data,
displacement_ref_felt,
local_displacement,
displacement_gradient);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment