Skip to content
GitLab
Commit e20b13f8 authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#820 SolidDeltaResidual: remove part that were related to bilinear contribution.

parent 0450bd10
Hide whitespace changes
Inline Side-by-side
Sources/ModelInstances/UnderDevelopment/Poromechanics/LocalVariationalOperatorInstances/SolidDeltaResidual.hpp
View file @ e20b13f8
......@@ -172,8 +172,7 @@ namespace HappyHeart
void ComputeAtQuadraturePoint(const InformationsAtQuadraturePoint& infos_at_quad_pt,
const LocalMatrix& tangent_matrix,
const LocalVector& rhs_part,
elementary_data_type& elementary_data,
OperatorNS::Nature what_must_be_built) const;
elementary_data_type& elementary_data) const;
//! Compute the derivates of W.
template<unsigned int DimensionT>
......@@ -194,7 +193,6 @@ namespace HappyHeart
const GeometricElt& geom_elt,
const Advanced::RefFEltInLocalOperator& ref_felt,
const std::vector<double>& local_displacement,
OperatorNS::Nature what_must_be_built,
const unsigned int mesh_dimension);
//! Helper class to manage computation of Cauchy-Green tensor.
......@@ -211,8 +209,7 @@ namespace HappyHeart
void PrepareInternalDataForQuadraturePointForDimension(const InformationsAtQuadraturePoint& infos_at_quad_pt,
const GeometricElt& geom_elt,
const Advanced::RefFEltInLocalOperator& ref_felt,
const std::vector<double>& local_displacement,
OperatorNS::Nature what_must_be_built);
const std::vector<double>& local_displacement);
private:
......
Sources/ModelInstances/UnderDevelopment/Poromechanics/LocalVariationalOperatorInstances/SolidDeltaResidual.hxx
View file @ e20b13f8
......@@ -124,9 +124,6 @@ namespace HappyHeart
const auto& local_displacement = GetFormerLocalDisplacement();
OperatorNS::Nature what_must_be_built = OperatorNS::Nature::nonlinear; // both matrix and vector must be built.
// see #533.
auto& tangent_matrix = this->matrix_parent::template GetLocalMatrix<EnumUnderlyingType(LocalMatrixIndex::tangent_matrix)>();
const auto& rhs_part = this->vector_parent::template GetLocalVector<EnumUnderlyingType(LocalVectorIndex::rhs_part)>();
......@@ -143,14 +140,12 @@ namespace HappyHeart
geom_elt,
elementary_data.GetRefFElt(this->GetNthUnknown()),
local_displacement,
what_must_be_built,
mesh_dimension);
ComputeAtQuadraturePoint(infos_at_quad_pt,
tangent_matrix,
rhs_part,
elementary_data,
what_must_be_built);
elementary_data);
}
}
......@@ -161,8 +156,7 @@ namespace HappyHeart
::ComputeAtQuadraturePoint(const InformationsAtQuadraturePoint& infos_at_quad_pt,
const LocalMatrix& tangent_matrix,
const LocalVector& rhs_part,
elementary_data_type& elementary_data,
OperatorNS::Nature what_must_be_built) const
elementary_data_type& elementary_data) const
{
const auto Ncomponent = elementary_data.GetGeomEltDimension();
......@@ -175,79 +169,27 @@ namespace HappyHeart
const auto& ref_felt = elementary_data.GetRefFElt(this->GetNthUnknown());
const int Nnode = static_cast<int>(elementary_data.Nnode());
if (what_must_be_built != OperatorNS::Nature::linear)
{
auto& gradient_based_block = this->matrix_parent::template GetLocalMatrix<EnumUnderlyingType(LocalMatrixIndex::gradient_based_block)>();
auto& transposed_dphi = this->matrix_parent::template GetLocalMatrix<EnumUnderlyingType(LocalMatrixIndex::transposed_dphi)>();
Wrappers::Seldon::Transpose(dPhi, transposed_dphi);
// Matrix related calculation.
auto& matrix_result = elementary_data.GetNonCstMatrixResult();
// LocalMatrix dPhi_mult_gradient_based_block(dPhi.GetM(), static_cast<int>(Ncomponent));
auto& dPhi_mult_gradient_based_block =
this->matrix_parent::template GetLocalMatrix<EnumUnderlyingType(LocalMatrixIndex::dPhi_mult_gradient_based_block)>();
auto& block_contribution =
this->matrix_parent::template GetLocalMatrix<EnumUnderlyingType(LocalMatrixIndex::block_contribution)>();
for (unsigned int row_component = 0u; row_component < Ncomponent; ++row_component)
{
const auto row_first_index = static_cast<int>(ref_felt.GetIndexFirstDofInElementaryData(row_component));
for (unsigned int col_component = 0u; col_component < Ncomponent; ++col_component)
{
const auto col_first_index = static_cast<int>(ref_felt.GetIndexFirstDofInElementaryData(col_component));
::HappyHeart::Private::ExtractGradientBasedBlock(tangent_matrix,
row_component,
col_component,
gradient_based_block);
Seldon::Mlt(1., dPhi,
gradient_based_block,
dPhi_mult_gradient_based_block);
Seldon::Mlt(weight_meas,
dPhi_mult_gradient_based_block,
transposed_dphi,
block_contribution);
for (int row_node = 0; row_node < Nnode; ++row_node)
{
for (int col_node = 0u; col_node < Nnode; ++col_node)
matrix_result(row_first_index + row_node, col_first_index + col_node)
+= block_contribution(row_node, col_node);
}
}
}
}
const auto int_Ncomponent = static_cast<int>(Ncomponent);
if (what_must_be_built != OperatorNS::Nature::bilinear)
auto& vector_result = elementary_data.GetNonCstVectorResult();
for (unsigned int row_component = 0u; row_component < Ncomponent; ++row_component)
{
// Vector related calculation.
auto& vector_result = elementary_data.GetNonCstVectorResult();
for (unsigned int row_component = 0u; row_component < Ncomponent; ++row_component)
{
const auto dof_first_index = static_cast<int>(ref_felt.GetIndexFirstDofInElementaryData(row_component));
const auto component_first_index = static_cast<int>(row_component * Ncomponent);
const auto dof_first_index = static_cast<int>(ref_felt.GetIndexFirstDofInElementaryData(row_component));
const auto component_first_index = static_cast<int>(row_component * Ncomponent);
// 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 (int row_node = 0; row_node < Nnode; ++row_node)
{
double value = 0.;
for (int col = 0; col < int_Ncomponent; ++col)
value += dPhi(row_node, col) * rhs_part(col + component_first_index);
vector_result(dof_first_index + row_node) += value * weight_meas;
}
}
}
// 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 (int row_node = 0; row_node < Nnode; ++row_node)
{
double value = 0.;
for (int col = 0; col < int_Ncomponent; ++col)
value += dPhi(row_node, col) * rhs_part(col + component_first_index);
vector_result(dof_first_index + row_node) += value * weight_meas;
}
}
}
......@@ -283,7 +225,6 @@ namespace HappyHeart
const GeometricElt& geom_elt,
const Advanced::RefFEltInLocalOperator& ref_felt,
const std::vector<double>& local_displacement,
OperatorNS::Nature what_must_be_built,
const unsigned int mesh_dimension)
{
switch(mesh_dimension)
......@@ -293,8 +234,7 @@ namespace HappyHeart
this->PrepareInternalDataForQuadraturePointForDimension<2>(infos_at_quad_pt,
geom_elt,
ref_felt,
local_displacement,
what_must_be_built);
local_displacement);
break;
}
case 3:
......@@ -302,8 +242,7 @@ namespace HappyHeart
this->PrepareInternalDataForQuadraturePointForDimension<3>(infos_at_quad_pt,
geom_elt,
ref_felt,
local_displacement,
what_must_be_built);
local_displacement);
break;
}
default:
......@@ -355,8 +294,7 @@ namespace HappyHeart
::PrepareInternalDataForQuadraturePointForDimension(const InformationsAtQuadraturePoint& infos_at_quad_pt,
const GeometricElt& geom_elt,
const Advanced::RefFEltInLocalOperator& ref_felt,
const std::vector<double>& local_displacement,
OperatorNS::Nature what_must_be_built)
const std::vector<double>& local_displacement)
{
auto& transposed_dphi = this->matrix_parent::template GetLocalMatrix<EnumUnderlyingType(LocalMatrixIndex::transposed_dphi)>();
Wrappers::Seldon::Transpose(infos_at_quad_pt.GetGradientFEltPhi(), transposed_dphi);
......@@ -402,14 +340,7 @@ namespace HappyHeart
dW,
d2W);
if (geom_elt.GetIndex() == 10000 || geom_elt.GetIndex() == 10000 || geom_elt.GetIndex() == 10000 || geom_elt.GetIndex() == 10000)
{
std::cout << "dW after hyper visco active" << std::endl;
for (int i = 0; i < dW.GetLength(); i++)
std::cout << std::scientific << std::setprecision(12) << dW(i) << "\t";
std::cout << std::endl;
}
# ifndef NDEBUG
Wrappers::Seldon::ThrowIfNanInside(d2W, __FILE__, __LINE__);
# endif // NDEBUG
......@@ -418,20 +349,8 @@ namespace HappyHeart
// Finally build the terms that are actually required.
// ===================================================================================
// Linear term.
switch(what_must_be_built)
{
case OperatorNS::Nature::linear:
case OperatorNS::Nature::nonlinear:
{
auto& rhs_part = this->vector_parent::template GetLocalVector<EnumUnderlyingType(LocalVectorIndex::rhs_part)>();
Seldon::Mlt(transposed_De, dW, rhs_part);
break;
}
case OperatorNS::Nature::bilinear:
break;
}
auto& rhs_part = this->vector_parent::template GetLocalVector<EnumUnderlyingType(LocalVectorIndex::rhs_part)>();
Seldon::Mlt(transposed_De, dW, rhs_part);
}
......
Supports Markdown
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