Commit 6766e26e authored by GILLES Sebastien's avatar GILLES Sebastien

#1498 - #1292 MixedSolidIncompressibility was not properly updated following...

#1498 - #1292 MixedSolidIncompressibility was not properly updated following transition from Seldon to Xtensor.
parent 6781ec98
......@@ -274,17 +274,12 @@ namespace MoReFEM::Advanced::LocalVariationalOperatorNS
auto& vector_result = elementary_data.GetNonCstVectorResult();
const unsigned int Ncomponent = displacement_ref_felt.GetMeshDimension();
const intNcomponent = static_cast<int>(Ncomponent);
const unsigned int Nnode_disp = displacement_ref_felt.Nnode();
const unsigned int Nnode_pressure = pressure_ref_felt.Nnode();
const int int_Nnode_disp = static_cast<int>(Nnode_disp);
const int int_Nnode_pressure = static_cast<int>(Nnode_pressure);
const unsigned int Nnode_test_disp = test_displacement_ref_felt.Nnode();
const unsigned int Nnode_test_pressure = test_pressure_ref_felt.Nnode();
const int int_Nnode_test_disp = static_cast<int>(Nnode_test_disp);
const int int_Nnode_test_pressure = static_cast<int>(Nnode_test_pressure);
const auto& quad_pt = infos_at_quad_pt.GetQuadraturePoint();
......@@ -361,19 +356,17 @@ namespace MoReFEM::Advanced::LocalVariationalOperatorNS
invariant_holder.Update(cauchy_green_tensor_at_quad_point, quad_pt, geom_elt);
double pressure_at_quad_point = 0.;
for (auto node_index = 0; node_index < int_Nnode_pressure; ++node_index)
for (auto node_index = 0ul; node_index < Nnode_pressure; ++node_index)
{
const unsigned int unsigned_node_index = static_cast<unsigned int>(node_index);
pressure_at_quad_point += pressure_phi(node_index) * local_pressure[unsigned_node_index];
pressure_at_quad_point += pressure_phi(node_index) * local_pressure[node_index];
}
// TODO : replace by template
const auto& dI3dC =
invariant_holder.GetFirstDerivativeWrtCauchyGreen(invariant_holder_type::invariants_first_derivative_index::dI3dC);
hydrostatic_stress.Copy(dI3dC);
diff_hydrostatic_stress_wrt_pres.Copy(hydrostatic_stress);
hydrostatic_stress = dI3dC;
diff_hydrostatic_stress_wrt_pres = hydrostatic_stress;
const double inv3 = Invariant3<FeltSpaceDimensionT>(cauchy_green_tensor_at_quad_point);
const double sqrt_inv3 = std::sqrt(inv3);
......@@ -396,11 +389,11 @@ namespace MoReFEM::Advanced::LocalVariationalOperatorNS
// Compute the new contribution to vector_result here.
// Product matrix vector is inlined here to avoid creation of an intermediate subset of \a rhs_part.
for (auto row_node = 0ul; row_node < int_Nnode_disp; ++row_node)
for (auto row_node = 0ul; row_node < Nnode_disp; ++row_node)
{
double value = 0.;
for (int col = 0; col <Ncomponent; ++col)
for (auto col = 0ul; col < Ncomponent; ++col)
value += dphi_test_displacement(row_node, col) * rhs_disp(col + component_first_index);
vector_result(dof_first_index + row_node) += value * weight_meas;
......@@ -409,7 +402,7 @@ namespace MoReFEM::Advanced::LocalVariationalOperatorNS
// Residual on pres
const auto dof_first_index_pres = static_cast<int>(Nnode_disp * Ncomponent);
for (auto row_node = 0ul; row_node < int_Nnode_pressure; ++row_node)
for (auto row_node = 0ul; row_node < Nnode_pressure; ++row_node)
{
vector_result(dof_first_index_pres + row_node) += rhs_pres * weight_meas * test_pressure_phi(row_node);
}
......@@ -453,7 +446,7 @@ namespace MoReFEM::Advanced::LocalVariationalOperatorNS
quad_pt,
geom_elt);
tangent_vector_pres_disp.Copy(tangent_vector_disp_pres);
tangent_vector_pres_disp = tangent_vector_disp_pres;
xt::noalias(tangent_vector_pres_disp) *= -2. * sqrt_inv3 * penalization_scd_deriv;
//std::cout << "tangent_vector_disp_pres"<< std::endl;
......@@ -482,9 +475,9 @@ namespace MoReFEM::Advanced::LocalVariationalOperatorNS
xt::noalias(block_contribution) = weight_meas * xt::linalg::dot(dphi_test_disp_mult_gradient_based_block,
transposed_dphi_displacement);
for (auto row_node = 0ul; row_node < int_Nnode_test_disp; ++row_node)
for (auto row_node = 0ul; row_node < Nnode_test_disp; ++row_node)
{
for (int col_node = 0u; col_node < int_Nnode_disp; ++col_node)
for (auto col_node = 0ul; col_node < Nnode_disp; ++col_node)
matrix_result(row_first_index + row_node, col_first_index + col_node)
+= block_contribution(row_node, col_node);
}
......@@ -508,18 +501,18 @@ namespace MoReFEM::Advanced::LocalVariationalOperatorNS
xt::noalias(dphi_test_disp_mult_column_matrix) = weight_meas * xt::linalg::dot(dphi_test_displacement,
column_matrix);
for (auto row_node = 0ul; row_node < int_Nnode_test_disp; ++row_node)
for (auto row_node = 0ul; row_node < Nnode_test_disp; ++row_node)
{
for (auto col_node = 0ul; col_node < int_Nnode_pressure; ++col_node)
for (auto col_node = 0ul; col_node < Nnode_pressure; ++col_node)
{
block_contribution_disp_pres(row_node, col_node) =
dphi_test_disp_mult_column_matrix(row_node, 0) * pressure_phi(col_node);
}
}
for (auto row_node = 0ul; row_node < int_Nnode_test_disp; ++row_node)
for (auto row_node = 0ul; row_node < Nnode_test_disp; ++row_node)
{
for (auto col_node = 0ul; col_node < int_Nnode_pressure; ++col_node)
for (auto col_node = 0ul; col_node < Nnode_pressure; ++col_node)
matrix_result(row_first_index_disp + row_node, col_first_index_pressure + col_node)
+= block_contribution_disp_pres(row_node, col_node);
}
......@@ -538,27 +531,27 @@ namespace MoReFEM::Advanced::LocalVariationalOperatorNS
xt::noalias(row_matrix_mult_transposed_dphi_disp) = weight_meas * xt::linalg::dot(row_matrix,
transposed_dphi_displacement);
for (auto row_node = 0ul; row_node < int_Nnode_test_pressure; ++row_node)
for (auto row_node = 0ul; row_node < Nnode_test_pressure; ++row_node)
{
for (auto col_node = 0ul; col_node < int_Nnode_disp; ++col_node)
for (auto col_node = 0ul; col_node < Nnode_disp; ++col_node)
{
block_contribution_pres_disp(row_node, col_node) = test_pressure_phi(row_node) *
row_matrix_mult_transposed_dphi_disp(0, col_node);
}
}
for (auto row_node = 0ul; row_node < int_Nnode_test_pressure; ++row_node)
for (auto row_node = 0ul; row_node < Nnode_test_pressure; ++row_node)
{
for (int col_node = 0u; col_node < int_Nnode_disp; ++col_node)
for (auto col_node = 0ul; col_node < Nnode_disp; ++col_node)
matrix_result(row_first_index_pressure + row_node, col_first_index_disp + col_node)
+= block_contribution_pres_disp(row_node, col_node);
}
}
// Tangent pres pres
for (auto row_node = 0ul; row_node < int_Nnode_test_pressure; ++row_node)
for (auto row_node = 0ul; row_node < Nnode_test_pressure; ++row_node)
{
for (int col_node = 0u; col_node < int_Nnode_pressure; ++col_node)
for (auto col_node = 0ul; col_node < Nnode_pressure; ++col_node)
matrix_result(row_first_index_pressure + row_node, col_first_index_pressure + col_node)
+= tangent_matrix_pres_pres * test_pressure_phi(row_node) * pressure_phi(col_node);
}
......
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