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

#820 Poromechanics: rewrite SolidOnFluidMesh with interpolator P1 -> P1b.

parent c0aa1ee0
Hide whitespace changes
Inline Side-by-side
Data/Lua/demo_input_poromechanics.lua
View file @ fa23bc60
......@@ -1038,36 +1038,6 @@ bulk_solid = {
} -- bulk_solid
-- Solid displacement on fluid
InitVertexMatchingInterpolator10 = {
-- Finite element space for which the dofs index will be associated to each vertex.
-- Expected format: VALUE
finite_element_space = 10,
-- Numbering subsetfor which the dofs index will be associated to each vertex.
-- Expected format: "VALUE"
numbering_subset = 14
} -- InitVertexMatchingInterpolator10
-- Solid displacement
InitVertexMatchingInterpolator20 = {
-- Finite element space for which the dofs index will be associated to each vertex.
-- Expected format: VALUE
finite_element_space = 20,
-- Numbering subsetfor which the dofs index will be associated to each vertex.
-- Expected format: "VALUE"
numbering_subset = 20
} -- InitVertexMatchingInterpolator20
interpolation_file = "${HOME}/Codes/HappyHeart/Data/Interpolation/Poromechanics/InterfaceFluidSolid.hhdata"
-- Initial value of porosity at each dof.
-- Expected format: VALUE
......
HappyHeart.xcodeproj/project.pbxproj
View file @ fa23bc60
......@@ -14442,7 +14442,7 @@
CLANG_WARN_DIRECT_OBJC_ISA_USAGE = YES_ERROR;
CLANG_WARN_OBJC_ROOT_CLASS = YES_ERROR;
CLANG_WARN_UNREACHABLE_CODE = YES;
CODE_SIGN_IDENTITY = "-";
CODE_SIGN_IDENTITY = "";
ENABLE_STRICT_OBJC_MSGSEND = YES;
GCC_C_LANGUAGE_STANDARD = c99;
GCC_NO_COMMON_BLOCKS = YES;
......@@ -14450,7 +14450,6 @@
"DEBUG=1",
"$(inherited)",
);
MTL_ENABLE_DEBUG_INFO = YES;
PRODUCT_NAME = "$(TARGET_NAME)";
};
name = Debug;
......@@ -14464,13 +14463,12 @@
CLANG_WARN_DIRECT_OBJC_ISA_USAGE = YES_ERROR;
CLANG_WARN_OBJC_ROOT_CLASS = YES_ERROR;
CLANG_WARN_UNREACHABLE_CODE = YES;
CODE_SIGN_IDENTITY = "-";
COPY_PHASE_STRIP = NO;
CODE_SIGN_IDENTITY = "";
COPY_PHASE_STRIP = YES;
ENABLE_NS_ASSERTIONS = NO;
ENABLE_STRICT_OBJC_MSGSEND = YES;
GCC_C_LANGUAGE_STANDARD = c99;
GCC_NO_COMMON_BLOCKS = YES;
MTL_ENABLE_DEBUG_INFO = NO;
PRODUCT_NAME = "$(TARGET_NAME)";
};
name = Release;
......@@ -14805,7 +14803,7 @@
COPY_PHASE_STRIP = NO;
DEAD_CODE_STRIPPING = YES;
DEBUG_INFORMATION_FORMAT = dwarf;
ENABLE_TESTABILITY = YES;
ENABLE_TESTABILITY = NO;
GCC_C_LANGUAGE_STANDARD = c99;
GCC_DYNAMIC_NO_PIC = NO;
GCC_ENABLE_OBJC_EXCEPTIONS = YES;
......@@ -14884,6 +14882,7 @@
PRESERVE_DEAD_CODE_INITS_AND_TERMS = NO;
SCAN_ALL_SOURCE_FILES_FOR_INCLUDES = YES;
SDKROOT = macosx;
STRIP_INSTALLED_PRODUCT = NO;
USER_HEADER_SEARCH_PATHS = Sources;
VALIDATE_PRODUCT = NO;
WARNING_CFLAGS = (
......
Sources/ModelInstances/UnderDevelopment/Poromechanics/ExplicitStepVariationalFormulation.cpp
View file @ fa23bc60
......@@ -148,7 +148,6 @@ namespace HappyHeart
modified_velocity_previous_iteration.Copy(velocity_previous_time_iteration, __FILE__, __LINE__);
// \todo #820 Should be an interpolation here!!!
Wrappers::Petsc::AXPY(-1.,
solid_on_fluid_mesh.GetHalfSumVelocity(),
modified_velocity_previous_iteration,
......
Sources/ModelInstances/UnderDevelopment/Poromechanics/InputParameterList.hpp
View file @ fa23bc60
......@@ -151,11 +151,6 @@ namespace HappyHeart
InputParameter::Fluid::Density,
InputParameter::Fluid::Viscosity,
InputParameter::InitVertexMatchingInterpolator<EnumUnderlyingType(InitVertexMatchingInterpolator::solid)>,
InputParameter::InitVertexMatchingInterpolator<EnumUnderlyingType(InitVertexMatchingInterpolator::fluid)>,
InputParameter::InterpolationFile,
// \todo #820 This parameter probably won't remain; keep it at the moment until it's completely sure...
InputParameter::PoromechanicsNS::Porosity,
......
Sources/ModelInstances/UnderDevelopment/Poromechanics/ModelInitialize.hxx
View file @ fa23bc60
......@@ -78,13 +78,19 @@ namespace HappyHeart
{
const auto& numbering_subset =
const auto& numbering_subset_p1 =
god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::solid_displacement_in_fluid_mesh));
const auto& numbering_subset_p1b =
god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity));
solid_on_fluid_mesh_ =
std::make_unique<Private::SolidOnFluidMesh>(input_parameter_data,
god_of_dof,
numbering_subset,
numbering_subset_p1,
numbering_subset_p1b,
GetSolidDisplacement(),
time_manager);
}
......
Sources/ModelInstances/UnderDevelopment/Poromechanics/Private/SolidOnFluidMesh.cpp
View file @ fa23bc60
......@@ -8,13 +8,15 @@
// Copyright © 2016 Inria. All rights reserved.
//
# include "ThirdParty/Wrappers/Petsc/Matrix/MatrixOperations.hpp"
#include "Core/NumberingSubset.hpp"
#include "Core/TimeManager/TimeManager.hpp"
#include "FiniteElement/FiniteElementSpace/GodOfDof.hpp"
#include "FiniteElement/FiniteElementSpace/Private/DofProgramWiseIndexListPerVertexCoordIndexManager.hpp"
#include "Operators/NonConformInterpolator/FromVertexMatching.hpp"
#include "Operators/ConformInterpolatorInstances/P1_to_P1b.hpp"
#include "ModelInstances/UnderDevelopment/Poromechanics/Private/SolidOnFluidMesh.hpp"
......@@ -31,72 +33,123 @@ namespace HappyHeart
namespace Private
{
namespace // anonymous
{
void ComputeInterpolationMatrix(const InputParameterList& input_parameter_data,
GlobalMatrix& out)
{
auto& init_vertex_matching_manager =
::HappyHeart::Private::DofProgramWiseIndexListPerVertexCoordIndexManager::GetInstance();
NonConformInterpolatorNS::FromVertexMatching vertex_matching
(input_parameter_data,
init_vertex_matching_manager.GetDofProgramWiseIndexListPerVertexCoordIndex(EnumUnderlyingType(InitVertexMatchingInterpolator::solid)),
init_vertex_matching_manager.GetDofProgramWiseIndexListPerVertexCoordIndex(EnumUnderlyingType(InitVertexMatchingInterpolator::fluid)));
}
} // namespace anonymous
SolidOnFluidMesh::SolidOnFluidMesh(const InputParameterList& input_parameter_data,
const GodOfDof& god_of_dof,
const NumberingSubset& numbering_subset,
const NumberingSubset& numbering_subset_p1,
const NumberingSubset& numbering_subset_p1b,
const GlobalVector& solid_displacement_on_solid_mesh,
const TimeManager& time_manager)
: time_manager_(time_manager),
solid_displacement_on_solid_mesh_(solid_displacement_on_solid_mesh)
{
former_displacement_ = std::make_unique<GlobalVector>(numbering_subset);
auto& former_disp = GetNonCstFormerDisplacement();
AllocateGlobalVector(god_of_dof, former_disp);
former_displacement_on_solid_mesh_ = std::make_unique<GlobalVector>(solid_displacement_on_solid_mesh);
former_displacement_on_solid_mesh_->ZeroEntries(__FILE__, __LINE__);
half_sum_velocity_on_solid_mesh_ = std::make_unique<GlobalVector>(solid_displacement_on_solid_mesh);
current_displacement_ = std::make_unique<GlobalVector>(numbering_subset_p1);
auto& current_displacement = GetNonCstCurrentDisplacement();
AllocateGlobalVector(god_of_dof, current_displacement);
half_sum_velocity_ = std::make_unique<GlobalVector>(numbering_subset_p1b);
auto& half_sum_velocity = GetNonCstHalfSumVelocity();
AllocateGlobalVector(god_of_dof, half_sum_velocity);
current_displacement_ = std::make_unique<GlobalVector>(former_disp);
half_sum_velocity_ = std::make_unique<GlobalVector>(former_disp);
InitInterpolationMatrix(input_parameter_data, god_of_dof);
}
void SolidOnFluidMesh::Update()
{
auto& former_disp = GetNonCstFormerDisplacement();
auto& current_disp = GetNonCstCurrentDisplacement();
// First update quantities on the solid mesh.
{
auto& former_disp = GetNonCstFormerDisplacementOnSolidMesh();
const auto& solid_displacement_on_solid_mesh = GetSolidDisplacementOnSolidMesh();
auto& half_sum_vel_on_solid_mesh = GetNonCstHalfSumVelocityOnSolidMesh();
// \todo #530 Use rather a swap when it's sorted out!
half_sum_vel_on_solid_mesh.Copy(former_disp, __FILE__, __LINE__);
Wrappers::Petsc::AXPY(1., solid_displacement_on_solid_mesh,
half_sum_vel_on_solid_mesh,
__FILE__, __LINE__);
const auto time_step = GetTimeManager().GetTimeStep();
assert(!NumericNS::IsZero(time_step));
half_sum_vel_on_solid_mesh.Scale(0.5, __FILE__, __LINE__);
// Update former displacement for next call.
// \todo #530 Use rather a swap when it's sorted out!
former_disp.Copy(solid_displacement_on_solid_mesh, __FILE__, __LINE__);
}
// \todo #530 Use rather a swap when it's sorted out!
former_disp.Copy(current_disp, __FILE__, __LINE__);
// \todo #820 Currently it's a mere copy, but later on it should rather be an interpolation!
current_disp.Copy(GetSolidDisplacementOnSolidMesh(), __FILE__, __LINE__);
// Then interpolate on the fluid mesh.
{
const auto& half_sum_vel_on_solid_mesh = GetHalfSumVelocityOnSolidMesh();
auto& half_sum_vel = GetNonCstHalfSumVelocity();
Wrappers::Petsc::MatMult(GetMatrixSolidToFluid(),
half_sum_vel_on_solid_mesh,
half_sum_vel,
__FILE__, __LINE__);
}
auto& half_sum_vel = GetNonCstHalfSumVelocity();
half_sum_vel.Copy(current_disp, __FILE__, __LINE__);
Wrappers::Petsc::AXPY(1., former_disp, half_sum_vel, __FILE__, __LINE__);
const auto time_step = GetTimeManager().GetTimeStep();
assert(!NumericNS::IsZero(time_step));
half_sum_vel.Scale(0.5 / time_step, __FILE__, __LINE__);
}
namespace // anonymous
void SolidOnFluidMesh::InitInterpolationMatrix(const InputParameterList& input_parameter_data,
const GodOfDof& god_of_dof)
{
// auto& init_vertex_matching_manager =
// ::HappyHeart::Private::DofProgramWiseIndexListPerVertexCoordIndexManager::GetInstance();
//
// NonConformInterpolatorNS::FromVertexMatching vertex_matching
// (input_parameter_data,
// init_vertex_matching_manager.GetDofProgramWiseIndexListPerVertexCoordIndex(EnumUnderlyingType(InitVertexMatchingInterpolator::solid)),
// init_vertex_matching_manager.GetDofProgramWiseIndexListPerVertexCoordIndex(EnumUnderlyingType(InitVertexMatchingInterpolator::fluid)));
//
decltype(auto) felt_space = god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::fluid));
} // namespace anonymous
decltype(auto) unknown_manager = UnknownManager::GetInstance();
decltype(auto) fluid_vel_unknown_ptr =
unknown_manager.GetUnknownPtr(EnumUnderlyingType(UnknownIndex::fluid_velocity));
decltype(auto) solid_disp_on_fluid_unknown_ptr =
unknown_manager.GetUnknownPtr(EnumUnderlyingType(UnknownIndex::solid_displacement));
decltype(auto) solid_disp_numbering_subset =
felt_space.GetNumberingSubset(*solid_disp_on_fluid_unknown_ptr);
decltype(auto) fluid_vel_numbering_subset =
felt_space.GetNumberingSubset(*fluid_vel_unknown_ptr);
matrix_solid_to_fluid_ = std::make_unique<GlobalMatrix>(fluid_vel_numbering_subset,
solid_disp_numbering_subset);
auto pairing = { std::make_pair(solid_disp_on_fluid_unknown_ptr, fluid_vel_unknown_ptr) };
ConformInterpolatorNS::P1_to_P1b interpolator_P1_to_P1b(felt_space,
solid_disp_numbering_subset,
felt_space,
fluid_vel_numbering_subset,
std::move(pairing));
interpolator_P1_to_P1b.Init();
GetNonCstMatrixSolidToFluid().CompleteCopy(interpolator_P1_to_P1b.GetInterpolationMatrix(),
__FILE__, __LINE__);
// Wrappers::Petsc::MatMatMult(
// vertex_matching.GetInterpolationMatrix(),
// GetNonCstMatrixSolidToFluid(),
// __FILE__, __LINE__);
}
} // namespace Private
......
Sources/ModelInstances/UnderDevelopment/Poromechanics/Private/SolidOnFluidMesh.hpp
View file @ fa23bc60
......@@ -52,6 +52,9 @@ namespace HappyHeart
* \brief Helper class which holds the data about solid displacement/velocity in fluid mesh.
*
* This is rather convenient to avoid passing too much arguments in variational formulation constructors.
*
* \attention This class currently assumes solid and fluid mesh are identical; should be modified when
* finer mesh will be used for fluid mesh (but true interpolator would be needed at that point in this case).
*/
class SolidOnFluidMesh
{
......@@ -75,13 +78,15 @@ namespace HappyHeart
* \param[in] input_parameter_data Objects which stores the content extracted from the input parameter
* file.
* \param[in] god_of_dof God of dof related to the fluid mesh.
* \param[in] numbering_subset Numbering subset upon which solid displacement is defined (on fluid mesh...)
* \param[in] numbering_subset_p1 Numbering subset upon which solid displacement is defined (on fluid god of dof...)
* \param[in] numbering_subset_p1b Numbering subset upon which fluid velocity is defined.
* \param[in] solid_displacement_on_solid_mesh Solid displacement as computed on the solid mesh.
* \param[in] time_manager Time manager.
*/
explicit SolidOnFluidMesh(const InputParameterList& input_parameter_data,
const GodOfDof& god_of_dof,
const NumberingSubset& numbering_subset,
const NumberingSubset& numbering_subset_p1,
const NumberingSubset& numbering_subset_p1b,
const GlobalVector& solid_displacement_on_solid_mesh,
const TimeManager& time_manager);
......@@ -119,15 +124,26 @@ namespace HappyHeart
//! Constant accessor to the current displacement.
const GlobalVector& GetCurrentDisplacement() const noexcept;
//! Constant accessor to the displacement of previous time iteration.
const GlobalVector& GetFormerDisplacement() const noexcept;
//! Constant accessor to the computation of vs{n - 1/2} (half sum of current and former displacement divided by time step).
const GlobalVector& GetHalfSumVelocity() const noexcept;
///@}
private:
/*!
* \brief Init the interpolation matrix that transforms solid displacement on solid mesh (P1) into
* P1b displacement on fluid mesh.
*
* \param[in] input_parameter_data Objects which stores the content extracted from the input parameter
* file.
* \param[in] god_of_dof God of dof related to the fluid mesh.
*/
void InitInterpolationMatrix(const InputParameterList& input_parameter_data,
const GodOfDof& god_of_dof);
private:
/// \name Accessors.
......@@ -148,6 +164,60 @@ namespace HappyHeart
//! Non constant accessor to the computation of vs{n - 1/2} (half sum of current and former displacement divided by time step).
GlobalVector& GetNonCstHalfSumVelocity() noexcept;
/*!
* \brief Constant accessor to the matrix that transforms a P1 displacement on solid mesh into a P1b
* displacement on fluid mesh.
*/
const GlobalMatrix& GetMatrixSolidToFluid() const noexcept;
/*!
* \brief Non constant accessor to the matrix that transforms a P1 displacement on solid mesh into a
* P1b displacement on fluid mesh.
*/
GlobalMatrix& GetNonCstMatrixSolidToFluid() noexcept;
/*!
* \brief Constant accessor to the an internal vector to store the half sum of solid displacement prior
* to the itnerpolation.
*/
const GlobalVector& GetHalfSumDisplacement() const noexcept;
/*!
* \brief Non constant accessor to the an internal vector to store the half sum of solid displacement
* prior to the itnerpolation.
*/
GlobalVector& GetNonCstHalfSumDisplacement() noexcept;
/*!
* \brief Constant accessor to the an internal vector to store the former solid displacement on the
* solid mesh.
*/
const GlobalVector& GetFormerDisplacementOnSolidMesh() const noexcept;
/*!
* \brief Non constant accessor to the an internal vector to store the former solid displacement on the
* solid mesh.
*/
GlobalVector& GetNonCstFormerDisplacementOnSolidMesh() noexcept;
/*!
* \brief Constant accessor to the an internal vector to store the half sum of former and current solid
* velocity on the solid mesh.
*/
const GlobalVector& GetHalfSumVelocityOnSolidMesh() const noexcept;
/*!
* \brief Non constant accessor to the an internal vector to store the half sum of former and current
* solid velocity on the solid mesh.
*/
GlobalVector& GetNonCstHalfSumVelocityOnSolidMesh() noexcept;
///@}
private:
......@@ -155,18 +225,23 @@ namespace HappyHeart
//! TimeManager.
const TimeManager& time_manager_;
//! Current displacement.
//! Current displacement on the fluid mesh.
GlobalVector::unique_ptr current_displacement_ = nullptr;
//! Displacement of previous time iteration.
GlobalVector::unique_ptr former_displacement_ = nullptr;
//! Computation of vs{n - 1/2} (half sum of current and former displacement divided by time step).
GlobalVector::unique_ptr half_sum_velocity_ = nullptr;
//! Solid displacement as computed on the solid mesh.
const GlobalVector& solid_displacement_on_solid_mesh_;
//! An internal vector to store the former solid displacement on the solid mesh.
GlobalVector::unique_ptr former_displacement_on_solid_mesh_ = nullptr;
//! An internal vector to store the half sum of former and current solid velocity on the solid mesh.
GlobalVector::unique_ptr half_sum_velocity_on_solid_mesh_ = nullptr;
//! Matrix that transforms a P1 displacement on solid mesh into a P1b displacement on fluid mesh.
GlobalMatrix::unique_ptr matrix_solid_to_fluid_ = nullptr;
};
......
Sources/ModelInstances/UnderDevelopment/Poromechanics/Private/SolidOnFluidMesh.hxx
View file @ fa23bc60
......@@ -36,20 +36,7 @@ namespace HappyHeart
return const_cast<GlobalVector&>(GetCurrentDisplacement());
}
inline const GlobalVector& SolidOnFluidMesh::GetFormerDisplacement() const noexcept
{
assert(!(!former_displacement_));
return *former_displacement_;
}
inline GlobalVector& SolidOnFluidMesh::GetNonCstFormerDisplacement() noexcept
{
return const_cast<GlobalVector&>(GetFormerDisplacement());
}
inline const GlobalVector& SolidOnFluidMesh::GetSolidDisplacementOnSolidMesh() const noexcept
{
return solid_displacement_on_solid_mesh_;
......@@ -75,6 +62,46 @@ namespace HappyHeart
}
inline const GlobalMatrix& SolidOnFluidMesh::GetMatrixSolidToFluid() const noexcept
{
assert(!(!matrix_solid_to_fluid_));
return *matrix_solid_to_fluid_;
}
inline GlobalMatrix& SolidOnFluidMesh::GetNonCstMatrixSolidToFluid() noexcept
{
return const_cast<GlobalMatrix&>(GetMatrixSolidToFluid());
}
inline const GlobalVector& SolidOnFluidMesh::GetFormerDisplacementOnSolidMesh() const noexcept
{
assert(!(!former_displacement_on_solid_mesh_));
return *former_displacement_on_solid_mesh_;
}
inline GlobalVector& SolidOnFluidMesh::GetNonCstFormerDisplacementOnSolidMesh() noexcept
{
return const_cast<GlobalVector&>(GetFormerDisplacementOnSolidMesh());
}
inline const GlobalVector& SolidOnFluidMesh::GetHalfSumVelocityOnSolidMesh() const noexcept
{
assert(!(!half_sum_velocity_on_solid_mesh_));
return *half_sum_velocity_on_solid_mesh_;
}
inline GlobalVector& SolidOnFluidMesh::GetNonCstHalfSumVelocityOnSolidMesh() noexcept
{
return const_cast<GlobalVector&>(GetHalfSumVelocityOnSolidMesh());
}
} // namespace Private
......
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