Commit 3cca95ce authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#820 Poromechanics: correctly initializes the fluid mass so that first...

#820 Poromechanics: correctly initializes the fluid mass so that first computed porosity matches Freefem's one.
parent 5f5f5f89
......@@ -129,9 +129,7 @@ namespace HappyHeart
const auto factor = quad_pt_factor * phi(node_index);
for (int component = 0; component < Ncomponent; ++component, dof_index += Nnode)
{
vector_result(dof_index) += factor;
}
}
}
}
......
......@@ -169,6 +169,9 @@ namespace HappyHeart
//! Accessor to mass vector.
const GlobalVector& GetFluidMassVector() const noexcept;
//! Non constant accessor to mass vector.
GlobalVector& GetNonCstFluidMassVector() noexcept;
//! Accessor to fluid density.
const ScalarParameter& GetFluidDensity() const noexcept;
......@@ -202,8 +205,12 @@ namespace HappyHeart
///@{
//! Initialize mass vector.
void InitializeMassVector();
/*!
* \brief Initialize mass vector.
*
* \param[in] phi_0 See namesake in Bruno's documentation.
*/
void InitializeMassVector(const double phi_0);
//! Initialize mass operator.
void InitializeMassOperator();
......
......@@ -34,9 +34,7 @@ namespace HappyHeart
template<class SolidVariationalFormulationPolicyT>
void Model<SolidVariationalFormulationPolicyT>::SupplFinalize()
{
// TODO: Put there what to do when all the time steps are done
}
{ }
template<class SolidVariationalFormulationPolicyT>
......
......@@ -72,6 +72,13 @@ namespace HappyHeart
}
template<class SolidVariationalFormulationPolicyT>
inline GlobalVector& Model<SolidVariationalFormulationPolicyT>::GetNonCstFluidMassVector() noexcept
{
return const_cast<GlobalVector&>(GetFluidMassVector());
}
template<class SolidVariationalFormulationPolicyT>
inline const ScalarParameter& Model<SolidVariationalFormulationPolicyT>::GetFluidDensity() const noexcept
{
......
......@@ -53,7 +53,6 @@ namespace HappyHeart
// Ordering is extremely important here!
InitializeMassVector();
InitParameters(input_parameter_data);
InitPorosity(input_parameter_data);
......@@ -191,6 +190,14 @@ namespace HappyHeart
throw Exception("Current model is restricted to a constant fluid density! (in fact a density defined "
"at dof would still be fine, but not at quadrature point...)",
__FILE__, __LINE__);
{
namespace ipl = Utilities::InputParameterListNS;
const double phi_0 = ipl::Extract<InputParameter::PoromechanicsNS::Porosity>::Value(input_parameter_data);
InitializeMassVector(phi_0);
}
const auto& felt_space = god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::fluid));
......@@ -278,12 +285,18 @@ namespace HappyHeart
template<class SolidVariationalFormulationPolicyT>
void Model<SolidVariationalFormulationPolicyT>::InitializeMassVector()
void Model<SolidVariationalFormulationPolicyT>::InitializeMassVector(const double phi_0)
{
decltype(auto) god_of_dof = parent::GetGodOfDof(EnumUnderlyingType(MeshIndex::fluid));
decltype(auto) numbering_subset = god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_scalar));
fluid_mass_vector_ = std::make_unique<GlobalVector>(numbering_subset);
AllocateGlobalVector(god_of_dof, *fluid_mass_vector_);
auto& fluid_mass_vector = GetNonCstFluidMassVector();
AllocateGlobalVector(god_of_dof, fluid_mass_vector);
fluid_mass_vector.SetUniformValue(GetFluidDensity().GetConstantValue() * phi_0, __FILE__, __LINE__);
}
......
......@@ -238,6 +238,14 @@ namespace HappyHeart
}
void Vector::SetUniformValue(PetscScalar value, const char* invoking_file, int invoking_line)
{
int error_code = VecSet(Internal(), value);
if (error_code)
throw ExceptionNS::Exception(error_code, "VecSetValue", invoking_file, invoking_line);
}
std::vector<PetscScalar> Vector::GetValues(const std::vector<PetscInt>& indexing,
const char* invoking_file, int invoking_line) const
......
......@@ -357,8 +357,6 @@ namespace HappyHeart
/*!
* \brief Add or modify one value inside a Petsc vector.
*
* VecAssembly is also called afterwards to prevent cache problems (see Petsc documentation about VecSetValues)
*
* \param[in] index Index (program-wise) to be modified.
* \param[in] value Value to set.
* \param[in] insertOrAppend Petsc ADD_VALUES or INSERT_VALUES (see Petsc documentation).
......@@ -368,6 +366,17 @@ namespace HappyHeart
void SetValue(PetscInt index, PetscScalar value, InsertMode insertOrAppend, const char* invoking_file, int invoking_line);
/*!
* \brief Set the same value to all entries of the vector.
*
* \param[in] value Value to set.
* \param[in] invoking_file File that invoked the function or class; usually __FILE__.
* \param[in] invoking_line File that invoked the function or class; usually __LINE__.
*/
void SetUniformValue(PetscScalar value, const char* invoking_file, int invoking_line);
/*!
* \brief A wrapper over VecCopy, which assumed target already gets the right layout.
......
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