Skip to content

GitLab

Une MAJ de sécurité est nécessaire sur notre version actuelle. Elle sera effectuée lundi 02/08 entre 12h30 et 13h. L'interruption de service devrait durer quelques minutes (probablement moins de 5 minutes).

Commit 50e0e0c1 authored by Gautier Bureau's avatar Gautier Bureau Committed by GILLES Sebastien
Browse files

#884 New formulation working in sequential. Parallel still not sure.

parent 3db6e306
Hide whitespace changes
Inline Side-by-side
Showing with 279 additions and 327 deletions
HappyHeart.xcodeproj/xcshareddata/xcschemes/CardiacMechanicsPrestressDisplacementVelocityFormulation-Mpi.xcscheme 0 → 100644
View file @ 50e0e0c1
<?xml version="1.0" encoding="UTF-8"?>
<Scheme
LastUpgradeVersion = "0730"
version = "1.3">
<BuildAction
parallelizeBuildables = "YES"
buildImplicitDependencies = "YES">
<BuildActionEntries>
<BuildActionEntry
buildForTesting = "YES"
buildForRunning = "YES"
buildForProfiling = "YES"
buildForArchiving = "YES"
buildForAnalyzing = "YES">
<BuildableReference
BuildableIdentifier = "primary"
BlueprintIdentifier = "1343EC271CC4D942009F854F"
BuildableName = "CardiacMechanicsPrestressDisplacementVelocityFormulation"
BlueprintName = "CardiacMechanicsPrestressDisplacementVelocityFormulation"
ReferencedContainer = "container:HappyHeart.xcodeproj">
</BuildableReference>
</BuildActionEntry>
</BuildActionEntries>
</BuildAction>
<TestAction
buildConfiguration = "Debug"
selectedDebuggerIdentifier = "Xcode.DebuggerFoundation.Debugger.LLDB"
selectedLauncherIdentifier = "Xcode.DebuggerFoundation.Launcher.LLDB"
shouldUseLaunchSchemeArgsEnv = "YES">
<Testables>
</Testables>
<MacroExpansion>
<BuildableReference
BuildableIdentifier = "primary"
BlueprintIdentifier = "1343EC271CC4D942009F854F"
BuildableName = "CardiacMechanicsPrestressDisplacementVelocityFormulation"
BlueprintName = "CardiacMechanicsPrestressDisplacementVelocityFormulation"
ReferencedContainer = "container:HappyHeart.xcodeproj">
</BuildableReference>
</MacroExpansion>
<AdditionalOptions>
</AdditionalOptions>
</TestAction>
<LaunchAction
buildConfiguration = "Debug"
selectedDebuggerIdentifier = "Xcode.DebuggerFoundation.Debugger.LLDB"
selectedLauncherIdentifier = "Xcode.DebuggerFoundation.Launcher.LLDB"
launchStyle = "0"
useCustomWorkingDirectory = "NO"
ignoresPersistentStateOnLaunch = "NO"
debugDocumentVersioning = "YES"
debugServiceExtension = "internal"
allowLocationSimulation = "YES">
<PathRunnable
runnableDebuggingMode = "0"
FilePath = "/Users/Shared/LibraryVersions/clang/Openmpi/bin/orterun">
</PathRunnable>
<MacroExpansion>
<BuildableReference
BuildableIdentifier = "primary"
BlueprintIdentifier = "1343EC271CC4D942009F854F"
BuildableName = "CardiacMechanicsPrestressDisplacementVelocityFormulation"
BlueprintName = "CardiacMechanicsPrestressDisplacementVelocityFormulation"
ReferencedContainer = "container:HappyHeart.xcodeproj">
</BuildableReference>
</MacroExpansion>
<CommandLineArguments>
<CommandLineArgument
argument = "-np 4"
isEnabled = "YES">
</CommandLineArgument>
<CommandLineArgument
argument = "$(BUILT_PRODUCTS_DIR)/CardiacMechanicsPrestressDisplacementVelocityFormulation"
isEnabled = "YES">
</CommandLineArgument>
<CommandLineArgument
argument = "--input_parameters ${HOME}/Codes/HappyHeart/Data/Lua/demo_input_cardiac_mechanics_prestress_displacement_velocity_formulation.lua"
isEnabled = "YES">
</CommandLineArgument>
</CommandLineArguments>
<AdditionalOptions>
</AdditionalOptions>
</LaunchAction>
<ProfileAction
buildConfiguration = "Release"
shouldUseLaunchSchemeArgsEnv = "YES"
savedToolIdentifier = ""
useCustomWorkingDirectory = "NO"
debugDocumentVersioning = "YES">
<BuildableProductRunnable
runnableDebuggingMode = "0">
<BuildableReference
BuildableIdentifier = "primary"
BlueprintIdentifier = "1343EC271CC4D942009F854F"
BuildableName = "CardiacMechanicsPrestressDisplacementVelocityFormulation"
BlueprintName = "CardiacMechanicsPrestressDisplacementVelocityFormulation"
ReferencedContainer = "container:HappyHeart.xcodeproj">
</BuildableReference>
</BuildableProductRunnable>
</ProfileAction>
<AnalyzeAction
buildConfiguration = "Debug">
</AnalyzeAction>
<ArchiveAction
buildConfiguration = "Release"
revealArchiveInOrganizer = "YES">
</ArchiveAction>
</Scheme>
HappyHeart.xcodeproj/xcshareddata/xcschemes/CardiacMechanicsPrestressDisplacementVelocityFormulation.xcscheme
View file @ 50e0e0c1
......@@ -42,7 +42,7 @@
</AdditionalOptions>
</TestAction>
<LaunchAction
buildConfiguration = "Debug"
buildConfiguration = "Release"
selectedDebuggerIdentifier = "Xcode.DebuggerFoundation.Debugger.LLDB"
selectedLauncherIdentifier = "Xcode.DebuggerFoundation.Launcher.LLDB"
launchStyle = "0"
......
Sources/ModelInstances/CardiacMechanicsPrestress/DisplacementFormulation/VariationalFormulation.hpp
View file @ 50e0e0c1
......@@ -87,9 +87,12 @@ namespace HappyHeart
HyperelasticLawNS::CiarletGeymonat
>;
//! Alias to the viscoelasticity derivate policy used.
using ViscoelasticityDerivate = LocalVariationalOperatorNS::SecondPiolaKirchhoffStressTensorNS::ViscoelasticityPolicyNS::DerivatesWithRespectTo;
//! Alias to the viscoelasticity policy used.
using ViscoelasticityPolicy =
GlobalVariationalOperatorNS::SecondPiolaKirchhoffStressTensorNS::ViscoelasticityPolicyNS::Viscoelasticity;
GlobalVariationalOperatorNS::SecondPiolaKirchhoffStressTensorNS::ViscoelasticityPolicyNS::Viscoelasticity<ViscoelasticityDerivate::displacement_and_velocity>;
//! Alias to the active stress policy used.
using ActiveStressPolicy =
......
Sources/ModelInstances/CardiacMechanicsPrestress/DisplacementVelocityFormulation/VariationalFormulation.cpp
View file @ 50e0e0c1
......@@ -23,9 +23,7 @@ namespace HappyHeart
namespace CardiacMechanicsPrestressDisplacementVelocityFormulationNS
{
int increment_file = 0;
VariationalFormulation::VariationalFormulation(const Wrappers::Mpi& mpi,
const NumberingSubset& displacement_numbering_subset,
......@@ -154,7 +152,6 @@ namespace HappyHeart
mpi.Barrier();
{
// \todo #583: Temporary to make code work; quite stupid nonetheless!
const auto degree_of_exactness = 2;
auto&& default_rule = DetermineDefaultQuadratureRule(degree_of_exactness);
GetNonCstQuadratureRulePerTopology().emplace_back(std::move(default_rule));
......@@ -183,9 +180,6 @@ namespace HappyHeart
parent::AllocateSystemMatrix(monolithic_displacement_velocity_numbering_subset, monolithic_displacement_velocity_numbering_subset);
parent::AllocateSystemVector(monolithic_displacement_velocity_numbering_subset);
//parent::AllocateSystemMatrix(displacement_numbering_subset, monolithic_displacement_velocity_numbering_subset);
//parent::AllocateSystemMatrix(velocity_numbering_subset, monolithic_displacement_velocity_numbering_subset);
parent::AllocateSystemVector(electrical_activation_numbering_subset);
const auto& displacement_system_matrix = GetSystemMatrix(displacement_numbering_subset, displacement_numbering_subset);
......@@ -196,16 +190,10 @@ namespace HappyHeart
const auto& monolithic_displacement_velocity_system_rhs = GetSystemRhs(monolithic_displacement_velocity_numbering_subset);
//const auto& displacement_to_displacement_monolithic_displacement_velocity_system_matrix = GetSystemMatrix(displacement_numbering_subset, monolithic_displacement_velocity_numbering_subset);
//const auto& velocity_to_displacement_monolithic_displacement_velocity_system_matrix = GetSystemMatrix(velocity_numbering_subset, monolithic_displacement_velocity_numbering_subset);
vector_stiffness_residual_ = std::make_unique<GlobalVector>(displacement_system_rhs);
vector_following_pressure_residual_ = std::make_unique<GlobalVector>(displacement_system_rhs);
vector_newmark_residual_ = std::make_unique<GlobalVector>(displacement_system_rhs);
vector_displacement_at_newton_iteration_ = std::make_unique<GlobalVector>(displacement_system_rhs);
vector_velocity_at_newton_iteration_ = std::make_unique<GlobalVector>(velocity_system_rhs);
......@@ -243,12 +231,16 @@ namespace HappyHeart
matrix_intermediate_displacement_to_monolithic_ = std::make_unique<GlobalMatrix>(displacement_numbering_subset, monolithic_displacement_velocity_numbering_subset);
matrix_intermediate_displacement_to_monolithic_2_ = std::make_unique<GlobalMatrix>(displacement_numbering_subset, monolithic_displacement_velocity_numbering_subset);
matrix_intermediate_velocity_to_monolithic_ = std::make_unique<GlobalMatrix>(velocity_numbering_subset, monolithic_displacement_velocity_numbering_subset);
matrix_intermediate_velocity_to_monolithic_2_ = std::make_unique<GlobalMatrix>(velocity_numbering_subset, monolithic_displacement_velocity_numbering_subset);
matrix_monolithic_global_ = std::make_unique<GlobalMatrix>(monolithic_displacement_velocity_numbering_subset, monolithic_displacement_velocity_numbering_subset);
matrix_monolithic_global_2_ = std::make_unique<GlobalMatrix>(monolithic_displacement_velocity_numbering_subset, monolithic_displacement_velocity_numbering_subset);
vector_monolithic_global_ = std::make_unique<GlobalVector>(monolithic_displacement_velocity_system_rhs);
}
}
......@@ -270,7 +262,6 @@ namespace HappyHeart
namespace IPL = Utilities::InputParameterListNS;
// \todo #583: Temporary to make code work; quite stupid nonetheless!
const auto& default_quadrature_rule_set = GetQuadratureRulePerTopologyType();
const auto& displacement_numbering_subset = GetDisplacementNumberingSubset();
......@@ -338,105 +329,10 @@ namespace HappyHeart
GetDisplacementLinearStiffnessOperator().Assemble(std::make_tuple(std::ref(matrix)));
}
GetMatrixDisplacementLinearStiffness().View(MpiHappyHeart(), "/Volumes/Data/gautier/HappyHeart/Computation/Cardiac/DispVel/K_disp.m", __FILE__, __LINE__, PETSC_VIEWER_ASCII_MATLAB);
{
GlobalMatrixWithCoefficient matrix(GetNonCstMatrixVelocityLinearStiffness(), 1.);
GetVelocityLinearStiffnessOperator().Assemble(std::make_tuple(std::ref(matrix)));
}
GetMatrixVelocityLinearStiffness().View(MpiHappyHeart(), "/Volumes/Data/gautier/HappyHeart/Computation/Cardiac/DispVel/K_vel.m", __FILE__, __LINE__, PETSC_VIEWER_ASCII_MATLAB);
// Get The list of dofs of the BC.
const auto& displacement_numbering_subset = GetDisplacementNumberingSubset();
const auto& dirichlet_bc_list = GetEssentialBoundaryConditionList();
std::vector<unsigned int> dofs_bc_disp;
for (const auto& dirichlet_bc_ptr : dirichlet_bc_list)
{
assert(!(!dirichlet_bc_ptr));
const auto& dirichlet_bc = *dirichlet_bc_ptr;
if (dirichlet_bc.IsNumberingSubset(displacement_numbering_subset))
{
const auto& dof_list = dirichlet_bc.GetDofList();
int unsigned size = static_cast<unsigned int>(dof_list.size());
for (unsigned int i = 0 ; i < size ; ++i)
dofs_bc_disp.push_back(dof_list[i]->GetProcessorWiseOrGhostIndex(displacement_numbering_subset));
}
}
std::ofstream outputFile ("/Volumes/Data/gautier/HappyHeart/Computation/Cardiac/DispVel/dofs_bc_disp_cardiac.dat");
if (outputFile.is_open())
{
for (unsigned int i = 0 ; i < dofs_bc_disp.size() ; ++i)
{
outputFile << dofs_bc_disp[i] << std::endl;
}
outputFile.close();
}
const auto& velocity_numbering_subset = GetVelocityNumberingSubset();
std::vector<unsigned int> dofs_bc_vel;
for (const auto& dirichlet_bc_ptr : dirichlet_bc_list)
{
assert(!(!dirichlet_bc_ptr));
const auto& dirichlet_bc = *dirichlet_bc_ptr;
if (dirichlet_bc.IsNumberingSubset(velocity_numbering_subset))
{
const auto& dof_list = dirichlet_bc.GetDofList();
int unsigned size = static_cast<unsigned int>(dof_list.size());
for (unsigned int i = 0 ; i < size ; ++i)
dofs_bc_vel.push_back(dof_list[i]->GetProcessorWiseOrGhostIndex(velocity_numbering_subset));
}
}
std::ofstream outputFile2 ("/Volumes/Data/gautier/HappyHeart/Computation/Cardiac/DispVel/dofs_bc_vel_cardiac.dat");
if (outputFile2.is_open())
{
for (unsigned int i = 0 ; i < dofs_bc_vel.size() ; ++i)
{
outputFile2 << dofs_bc_vel[i] << std::endl;
}
outputFile2.close();
}
const auto& monolithic_numbering_subset = GetMonolithicDisplacementVelocityNumberingSubset();
std::vector<unsigned int> dofs_bc_disp_vel;
for (const auto& dirichlet_bc_ptr : dirichlet_bc_list)
{
assert(!(!dirichlet_bc_ptr));
const auto& dirichlet_bc = *dirichlet_bc_ptr;
const auto& dof_list = dirichlet_bc.GetDofList();
int unsigned size = static_cast<unsigned int>(dof_list.size());
for (unsigned int i = 0 ; i < size ; ++i)
dofs_bc_disp_vel.push_back(dof_list[i]->GetProcessorWiseOrGhostIndex(monolithic_numbering_subset));
}
std::ofstream outputFile3 ("/Volumes/Data/gautier/HappyHeart/Computation/Cardiac/DispVel/dofs_bc_disp_vel_cardiac.dat");
if (outputFile3.is_open())
{
for (unsigned int i = 0 ; i < dofs_bc_disp_vel.size() ; ++i)
{
outputFile3 << dofs_bc_disp_vel[i] << std::endl;
}
outputFile3.close();
}
}
void VariationalFormulation::ComputeApexElevation()
......@@ -587,14 +483,6 @@ namespace HappyHeart
monolithic_displacement_velocity_numbering_subset);
velocity_to_monolithic_displacement_velocity_interpolator_->Init();
GetDisplacementToMonolithicDisplacementVelocityInterpolator().GetInterpolationMatrix().View(MpiHappyHeart(), "/Volumes/Data/gautier/HappyHeart/Computation/Cardiac/DispVel/matrix_interpolation_disp_to_mono.m", __FILE__, __LINE__, PETSC_VIEWER_ASCII_MATLAB);
GetVelocityToMonolithicDisplacementVelocityInterpolator().GetInterpolationMatrix().View(MpiHappyHeart(), "/Volumes/Data/gautier/HappyHeart/Computation/Cardiac/DispVel/matrix_interpolation_vel_to_mono.m", __FILE__, __LINE__, PETSC_VIEWER_ASCII_MATLAB);
GetMonolithicDisplacementVelocityToDisplacementInterpolator().GetInterpolationMatrix().View(MpiHappyHeart(), "/Volumes/Data/gautier/HappyHeart/Computation/Cardiac/DispVel/matrix_interpolation_mono_to_disp.m", __FILE__, __LINE__, PETSC_VIEWER_ASCII_MATLAB);
GetMonolithicDisplacementVelocityToVelocityInterpolator().GetInterpolationMatrix().View(MpiHappyHeart(), "/Volumes/Data/gautier/HappyHeart/Computation/Cardiac/DispVel/matrix_interpolation_mono_to_vel.m", __FILE__, __LINE__, PETSC_VIEWER_ASCII_MATLAB);
}
......@@ -609,8 +497,6 @@ namespace HappyHeart
GetMassOperator().Assemble(std::make_tuple(std::ref(matrix)));
}
GetMatrixMass().View(MpiHappyHeart(), "/Volumes/Data/gautier/HappyHeart/Computation/Cardiac/DispVel/Mass.m", __FILE__, __LINE__, PETSC_VIEWER_ASCII_MATLAB);
{
const auto do_reuse_matrix = Wrappers::Petsc::DoReuseMatrix::no;
......@@ -634,33 +520,31 @@ namespace HappyHeart
__FILE__, __LINE__,
do_reuse_matrix);
}
GetMatrixMonolithicLinearStiffnessDisplacementDisplacement().View(MpiHappyHeart(), "/Volumes/Data/gautier/HappyHeart/Computation/Cardiac/DispVel/K_disp_disp_mono.m", __FILE__, __LINE__, PETSC_VIEWER_ASCII_MATLAB);
{
const auto do_reuse_matrix = Wrappers::Petsc::DoReuseMatrix::no;
const auto& matrix_velocity_linear_stiffness = GetMatrixVelocityLinearStiffness();
auto& matrix_intermediate_displacement_to_monolithic = GetNonCstMatrixIntermediateDisplacementToMonolithic();
auto& matrix_intermediate_velocity_to_monolithic = GetNonCstMatrixIntermediateVelocityToMonolithic();
auto& matrix_monolithic_linear_stiffness_displacement_velocity = GetNonCstMatrixMonolithicLinearStiffnessDisplacementVelocity();
const auto& monolithic_displacement_velocity_to_velocity_interpolator = GetMonolithicDisplacementVelocityToVelocityInterpolator();
const auto& monolithic_displacement_velocity_to_displacement_interpolator = GetMonolithicDisplacementVelocityToDisplacementInterpolator();
// Yes because the structure of the matrix is allocated.
// No because the structure of the matrix is not yet allocated.
Wrappers::Petsc::MatMatMult(matrix_velocity_linear_stiffness,
monolithic_displacement_velocity_to_velocity_interpolator.GetInterpolationMatrix(),
matrix_intermediate_displacement_to_monolithic,
matrix_intermediate_velocity_to_monolithic,
__FILE__, __LINE__,
Wrappers::Petsc::DoReuseMatrix::yes);
do_reuse_matrix);
// No because the structure of the matrix is not yet allocated.
Wrappers::Petsc::MatTransposeMatMult(monolithic_displacement_velocity_to_displacement_interpolator.GetInterpolationMatrix(),
matrix_intermediate_displacement_to_monolithic,
matrix_intermediate_velocity_to_monolithic,
matrix_monolithic_linear_stiffness_displacement_velocity,
__FILE__, __LINE__,
Wrappers::Petsc::DoReuseMatrix::no);
do_reuse_matrix);
}
GetMatrixMonolithicLinearStiffnessDisplacementVelocity().View(MpiHappyHeart(), "/Volumes/Data/gautier/HappyHeart/Computation/Cardiac/DispVel/K_disp_vel_mono.m", __FILE__, __LINE__, PETSC_VIEWER_ASCII_MATLAB);
}
......@@ -737,28 +621,15 @@ namespace HappyHeart
Wrappers::Petsc::AXPY(1., GetVectorVelocityAtNewtonIteration(), midpoint_velocity, __FILE__, __LINE__);
midpoint_velocity.Scale(0.5, __FILE__, __LINE__);
midpoint_velocity.UpdateGhosts(__FILE__, __LINE__);
// Midpoint disp in vel NS
const auto& const_midpoint_monolithic_displacement_velocity = GetVectorMidpointMonolithicDisplacementVelocity();
const auto& monolithic_displacement_velocity_to_displacement_interpolator = GetMonolithicDisplacementVelocityToDisplacementInterpolator();
auto& midpoint_position_in_velocity_numbering_subset = GetNonCstVectorMidpointPositionInVelocityNumberingSubset();
Wrappers::Petsc::MatMult(monolithic_displacement_velocity_to_displacement_interpolator.GetInterpolationMatrix(),
const_midpoint_monolithic_displacement_velocity,
midpoint_position_in_velocity_numbering_subset,
__FILE__, __LINE__);
midpoint_position_in_velocity_numbering_subset.Copy(GetVectorMidpointPosition(), __FILE__, __LINE__);
midpoint_position_in_velocity_numbering_subset.UpdateGhosts(__FILE__, __LINE__);
// Midpoint vel in disp NS
const auto& monolithic_displacement_velocity_to_velocity_interpolator = GetMonolithicDisplacementVelocityToVelocityInterpolator();
auto& midpoint_velocity_in_displacement_numbering_subset = GetNonCstVectorMidpointVelocityInDisplacementNumberingSubset();
Wrappers::Petsc::MatMult(monolithic_displacement_velocity_to_velocity_interpolator.GetInterpolationMatrix(),
const_midpoint_monolithic_displacement_velocity,
midpoint_velocity_in_displacement_numbering_subset,
__FILE__, __LINE__);
midpoint_velocity_in_displacement_numbering_subset.Copy(GetVectorMidpointVelocity(), __FILE__, __LINE__);
midpoint_velocity_in_displacement_numbering_subset.UpdateGhosts(__FILE__, __LINE__);
if (newton_iteration == 0)
......@@ -766,17 +637,6 @@ namespace HappyHeart
GetNonCstVectorDisplacementPreviousTimeIteration().UpdateGhosts(__FILE__, __LINE__);
GetNonCstVectorVelocityPreviousTimeIteration().UpdateGhosts(__FILE__, __LINE__);
// Midpoint disp in vel NS
//auto& midpoint_position_in_velocity_numbering_subset = GetNonCstVectorMidpointPositionInVelocityNumberingSubset();
//midpoint_position_in_velocity_numbering_subset.Copy(GetVectorMidpointPosition(), __FILE__, __LINE__);
//midpoint_position_in_velocity_numbering_subset.UpdateGhosts(__FILE__, __LINE__);
// Midpoint vel in disp NS
//auto& midpoint_velocity_in_displacement_numbering_subset = GetNonCstVectorMidpointVelocityInDisplacementNumberingSubset();
//midpoint_velocity_in_displacement_numbering_subset.Copy(GetVectorMidpointVelocity(), __FILE__, __LINE__);
//midpoint_velocity_in_displacement_numbering_subset.UpdateGhosts(__FILE__, __LINE__);
}
......@@ -965,6 +825,8 @@ namespace HappyHeart
GetViscoelasticityOperator().Assemble(std::make_tuple(std::ref(mat), std::ref(vec)),
displacement_vector);
velocity_rhs.ZeroEntries(__FILE__, __LINE__);
}
}
......@@ -1005,7 +867,7 @@ namespace HappyHeart
velocity_rhs.ZeroEntries(__FILE__, __LINE__);
monolithic_rhs.ZeroEntries(__FILE__, __LINE__);
// Complete residual on vel R of (Rd R);
// Complete residual on velocity R of (Rd R);
const auto& stiffness_residual = GetVectorStiffnessResidual();
Wrappers::Petsc::AXPY(1.,
......@@ -1016,7 +878,7 @@ namespace HappyHeart
const auto& time_manager = GetTimeManager();
const double time_step = time_manager.GetTimeStep();
const auto& velocity_previous_time_iteration = GetVectorVelocityPreviousTimeIteration();
auto& diff_displacement = GetNonCstVectorDiffDisplacement(); // In reality here it is the acceleration
auto& diff_displacement = GetNonCstVectorDiffDisplacement(); // In reality here it is the acceleration.
diff_displacement.Copy(GetVectorVelocityAtNewtonIteration(), __FILE__, __LINE__);
Wrappers::Petsc::AXPY(-1.,
......@@ -1031,7 +893,7 @@ namespace HappyHeart
velocity_rhs,
velocity_rhs,
__FILE__, __LINE__);
const auto& velocity_to_monolithic_displacement_velocity_interpolator = GetVelocityToMonolithicDisplacementVelocityInterpolator();
Wrappers::Petsc::MatMult(velocity_to_monolithic_displacement_velocity_interpolator.GetInterpolationMatrix(),
......@@ -1039,9 +901,7 @@ namespace HappyHeart
monolithic_rhs,
__FILE__, __LINE__);
GetSystemRhs(velocity_numbering_subset).Print<MpiScale::processor_wise>(MpiHappyHeart(), "/Volumes/Data/gautier/HappyHeart/Computation/Cardiac/DispVel/R/R_" + std::to_string(EnumUnderlyingType(GetTimeManager().GetStaticOrDynamic())) + "_" + std::to_string(increment_file) + "_" + std::to_string(GetSnes().GetSnesIteration(__FILE__, __LINE__)) + ".hhdata", __FILE__, __LINE__);
// Complete residual on vel Rd of (Rd R);
// Complete residual on displacement Rd of (Rd R);
const auto& displacement_previous_time_iteration = GetVectorDisplacementPreviousTimeIteration();
const auto& displacement_at_newton_iteration = GetVectorDisplacementAtNewtonIteration();
const auto& midpoint_velocity = GetVectorMidpointVelocity();
......@@ -1068,10 +928,9 @@ namespace HappyHeart
displacement_rhs,
__FILE__, __LINE__);
GetSystemRhs(displacement_numbering_subset).Print<MpiScale::processor_wise>(MpiHappyHeart(), "/Volumes/Data/gautier/HappyHeart/Computation/Cardiac/DispVel/Rd/Rd_" + std::to_string(EnumUnderlyingType(GetTimeManager().GetStaticOrDynamic())) + "_" + std::to_string(increment_file) + "_" + std::to_string(GetSnes().GetSnesIteration(__FILE__, __LINE__)) + ".hhdata", __FILE__, __LINE__);
const auto& displacement_to_monolithic_displacement_velocity_interpolator = GetDisplacementToMonolithicDisplacementVelocityInterpolator();
auto& vector_monolithic_global = GetNonCstVectorMonolithicGlobal();
vector_monolithic_global.ZeroEntries(__FILE__, __LINE__);
Wrappers::Petsc::MatMult(displacement_to_monolithic_displacement_velocity_interpolator.GetInterpolationMatrix(),
displacement_rhs,
......@@ -1082,16 +941,17 @@ namespace HappyHeart
GetVectorMonolithicGlobal(),
monolithic_rhs,
__FILE__, __LINE__);
GetSystemRhs(monolithic_displacement_velocity_numbering_subset).Print<MpiScale::processor_wise>(MpiHappyHeart(), "/Volumes/Data/gautier/HappyHeart/Computation/Cardiac/DispVel/R_gen/R_gen_" + std::to_string(EnumUnderlyingType(GetTimeManager().GetStaticOrDynamic())) + "_" + std::to_string(increment_file) + "_" + std::to_string(GetSnes().GetSnesIteration(__FILE__, __LINE__)) + ".hhdata", __FILE__, __LINE__);
ApplyEssentialBoundaryCondition<OperatorNS::Nature::linear>(monolithic_displacement_velocity_numbering_subset,
monolithic_displacement_velocity_numbering_subset);
ApplyEssentialBoundaryCondition<OperatorNS::Nature::linear>(displacement_numbering_subset,
displacement_numbering_subset);
ApplyEssentialBoundaryCondition<OperatorNS::Nature::linear>(velocity_numbering_subset,
velocity_numbering_subset);
TangentCheck(monolithic_rhs);
TangentCheck(velocity_rhs);
}
void VariationalFormulation::ComputeStaticTangent()
......@@ -1116,28 +976,23 @@ namespace HappyHeart
void VariationalFormulation::ComputeDynamicTangent()
{
const auto do_reuse_matrix_yes = Wrappers::Petsc::DoReuseMatrix::yes;
{
static bool is_first_call = true;
const auto do_reuse_matrix = is_first_call
? Wrappers::Petsc::DoReuseMatrix::no
: Wrappers::Petsc::DoReuseMatrix::yes;
const auto& displacement_numbering_subset = GetDisplacementNumberingSubset();
const auto& velocity_numbering_subset = GetVelocityNumberingSubset();
const auto& monolithic_displacement_velocity_numbering_subset = GetMonolithicDisplacementVelocityNumberingSubset();
auto& displacement_system_matrix = GetNonCstSystemMatrix(displacement_numbering_subset, displacement_numbering_subset);
auto& velocity_system_matrix = GetNonCstSystemMatrix(velocity_numbering_subset, velocity_numbering_subset);
auto& monolithic_displacement_velocity_system_matrix = GetNonCstSystemMatrix(monolithic_displacement_velocity_numbering_subset, monolithic_displacement_velocity_numbering_subset);
auto& intermediate_displacement_to_monolithic = GetNonCstMatrixIntermediateDisplacementToMonolithic();
auto& intermediate_velocity_to_monolithic = GetNonCstMatrixIntermediateVelocityToMonolithic();
auto& intermediate_displacement_to_monolithic_2 = GetNonCstMatrixIntermediateDisplacementToMonolithic2();
auto& intermediate_velocity_to_monolithic_2 = GetNonCstMatrixIntermediateVelocityToMonolithic2();
auto& matrix_monolithic_global = GetNonCstMatrixMonolithicGlobal();
auto& matrix_monolithic_global_2 = GetNonCstMatrixMonolithicGlobal2();
//displacement_system_matrix.ZeroEntries(__FILE__, __LINE__);
velocity_system_matrix.ZeroEntries(__FILE__, __LINE__);
monolithic_displacement_velocity_system_matrix.ZeroEntries(__FILE__, __LINE__);
......@@ -1148,10 +1003,7 @@ namespace HappyHeart
const auto& time_manager = GetTimeManager();
const double time_step = time_manager.GetTimeStep();
displacement_system_matrix.Copy(GetMatrixTangentStiffness(), __FILE__, __LINE__);
displacement_system_matrix.Scale(0.5, __FILE__, __LINE__);
displacement_system_matrix.View(MpiHappyHeart(), "/Volumes/Data/gautier/HappyHeart/Computation/Cardiac/DispVel/dK_dY12/dK_dY12_" + std::to_string(EnumUnderlyingType(GetTimeManager().GetStaticOrDynamic())) + "_" + std::to_string(increment_file) + "_" + std::to_string(GetSnes().GetSnesIteration(__FILE__, __LINE__)) + ".m", __FILE__, __LINE__, PETSC_VIEWER_ASCII_MATLAB);
GetNonCstMatrixTangentStiffness().Scale(0.5, __FILE__, __LINE__);
velocity_system_matrix.Copy(GetMatrixMass(), __FILE__, __LINE__);
velocity_system_matrix.Scale(1. / time_step, __FILE__, __LINE__);
......@@ -1161,83 +1013,61 @@ namespace HappyHeart
velocity_system_matrix,
__FILE__, __LINE__);
GetMatrixTangentViscoelasticity().View(MpiHappyHeart(), "/Volumes/Data/gautier/HappyHeart/Computation/Cardiac/DispVel/dC_dYp12/dC_dYp12_" + std::to_string(EnumUnderlyingType(GetTimeManager().GetStaticOrDynamic())) + "_" + std::to_string(increment_file) + "_" + std::to_string(GetSnes().GetSnesIteration(__FILE__, __LINE__)) + ".m", __FILE__, __LINE__, PETSC_VIEWER_ASCII_MATLAB);
// Add in the global matrix: 1/2*dK_dY12.
Wr