Skip to content
GitLab
Menu
Projects
Groups
Snippets
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
MoReFEM
CoreLibrary
MoReFEM
Commits
e4f12a02
Commit
e4f12a02
authored
Oct 27, 2016
by
GILLES Sebastien
Browse files
#1022
Make ExplicitStep a template class and add new fluid pressure data to it.
parent
518301ad
Changes
11
Expand all
Hide whitespace changes
Inline
Side-by-side
HappyHeart.xcodeproj/project.pbxproj
View file @
e4f12a02
...
...
@@ -451,7 +451,6 @@
BE17E7C01D7EA9E300A1AE6A /* FluidVelocity.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BE91D7291D67248C00EFD107 /* FluidVelocity.cpp */; };
BE17E7C11D7EA9E300A1AE6A /* Porosity.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BE60DBFF1C7F4222007B334C /* Porosity.cpp */; };
BE17E7C31D7EA9E300A1AE6A /* main_st_venant_kirchhoff.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BE0B0BD51C6A19130038A9B9 /* main_st_venant_kirchhoff.cpp */; };
BE17E7C41D7EA9E300A1AE6A /* VariationalFormulation.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BEA78B911CE0925D00A6A185 /* VariationalFormulation.cpp */; };
BE17E7C51D7EA9E300A1AE6A /* SolidVelocity.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BE8B86C91D6AEFF100F48E88 /* SolidVelocity.cpp */; };
BE17E7C61D7EA9E300A1AE6A /* GradOnGradientBasedElasticityTensor.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BE669E171C7C9FCE00C521CA /* GradOnGradientBasedElasticityTensor.cpp */; };
BE17E7C71D7EA9E300A1AE6A /* LoadOnSolid.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BEDAF9CF1D17F51600D5AA0D /* LoadOnSolid.cpp */; };
...
...
@@ -1190,7 +1189,6 @@
BE6FB2111D8C3BD800B9F0FB /* Porosity.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BE8B86CD1D6B47B400F48E88 /* Porosity.cpp */; };
BE6FB2121D8C3BD800B9F0FB /* FluidPressure.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BE15E40E1D6DC67E0046E1B2 /* FluidPressure.cpp */; };
BE6FB2131D8C3BF300B9F0FB /* VariationalFormulation.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BEA78BB31CE092CA00A6A185 /* VariationalFormulation.cpp */; };
BE6FB2141D8C3BF300B9F0FB /* VariationalFormulation.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BEA78B911CE0925D00A6A185 /* VariationalFormulation.cpp */; };
BE6FB2201D8C3BF300B9F0FB /* Ale.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BEF183251C6DE3A3008A6F1E /* Ale.cpp */; };
BE6FB2211D8C3BF300B9F0FB /* GradOnGradientBasedElasticityTensor.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BE669E171C7C9FCE00C521CA /* GradOnGradientBasedElasticityTensor.cpp */; };
BE6FB2221D8C3BF300B9F0FB /* Porosity.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BE60DBFB1C7F417C007B334C /* Porosity.cpp */; };
...
...
@@ -6237,7 +6235,6 @@
BEA355D717D0971500FB643B /* OpsFunction.hxx */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = OpsFunction.hxx; sourceTree = "<group>"; };
BEA4FC4418214D8F002B2EA1 /* VariationalFormulation.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; lineEnding = 0; path = VariationalFormulation.hpp; sourceTree = "<group>"; xcLanguageSpecificationIdentifier = xcode.lang.cpp; };
BEA4FC4518214D8F002B2EA1 /* VariationalFormulation.hxx */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; lineEnding = 0; path = VariationalFormulation.hxx; sourceTree = "<group>"; xcLanguageSpecificationIdentifier = xcode.lang.cpp; };
BEA78B911CE0925D00A6A185 /* VariationalFormulation.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = VariationalFormulation.cpp; sourceTree = "<group>"; };
BEA78B921CE0925D00A6A185 /* VariationalFormulation.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = VariationalFormulation.hpp; sourceTree = "<group>"; };
BEA78B931CE0925D00A6A185 /* VariationalFormulation.hxx */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = VariationalFormulation.hxx; sourceTree = "<group>"; };
BEA78BB31CE092CA00A6A185 /* VariationalFormulation.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = VariationalFormulation.cpp; sourceTree = "<group>"; };
...
...
@@ -11242,7 +11239,6 @@
BEA78B901CE0925D00A6A185 /* ExplicitStep */ = {
isa = PBXGroup;
children = (
BEA78B911CE0925D00A6A185 /* VariationalFormulation.cpp */,
BEA78B921CE0925D00A6A185 /* VariationalFormulation.hpp */,
BEA78B931CE0925D00A6A185 /* VariationalFormulation.hxx */,
);
...
...
@@ -15473,7 +15469,6 @@
BE17E7C01D7EA9E300A1AE6A /* FluidVelocity.cpp in Sources */,
BE17E7C11D7EA9E300A1AE6A /* Porosity.cpp in Sources */,
BE17E7C31D7EA9E300A1AE6A /* main_st_venant_kirchhoff.cpp in Sources */,
BE17E7C41D7EA9E300A1AE6A /* VariationalFormulation.cpp in Sources */,
BE17E7C51D7EA9E300A1AE6A /* SolidVelocity.cpp in Sources */,
BE71B9371DB75AE5005DCB59 /* T33.cpp in Sources */,
BE17E7C61D7EA9E300A1AE6A /* GradOnGradientBasedElasticityTensor.cpp in Sources */,
...
...
@@ -15823,7 +15818,6 @@
BE71B9361DB75AE5005DCB59 /* T33.cpp in Sources */,
BE6FB22F1D8C3BF300B9F0FB /* PenalizationPorosity.cpp in Sources */,
BE6FB22A1D8C3BF300B9F0FB /* FluidSourceTerm.cpp in Sources */,
BE6FB2141D8C3BF300B9F0FB /* VariationalFormulation.cpp in Sources */,
BE71B92A1DB75AE5005DCB59 /* T33.cpp in Sources */,
BE6FB20D1D8C3BD800B9F0FB /* SolidDisplacement.cpp in Sources */,
BE6FB2261D8C3BF300B9F0FB /* GradOnGradientBasedElasticityTensor.cpp in Sources */,
Sources/ModelInstances/UnderDevelopment/Poromechanics/Data/FluidVelocity.cpp
View file @
e4f12a02
...
...
@@ -100,20 +100,7 @@ namespace HappyHeart
}
void
FluidVelocity
::
PrepareNextTimeIteration
(
const
ExplicitStepVariationalFormulation
&
explicit_varf
)
{
auto
&
current
=
GetNonCstCurrent
();
// \todo #820 Swap...
current
.
Copy
(
GetNew
(),
__FILE__
,
__LINE__
);
Wrappers
::
Petsc
::
MatMult
(
explicit_varf
.
GetMassMatrixWithPorosity
(),
current
,
GetNonCstMassTimesCurrent
(),
__FILE__
,
__LINE__
);
}
}
// namespace DataNS
...
...
Sources/ModelInstances/UnderDevelopment/Poromechanics/Data/FluidVelocity.hpp
View file @
e4f12a02
...
...
@@ -37,6 +37,7 @@ namespace HappyHeart
// ============================
template
<
class
HyperelasticLawT
>
class
ExplicitStepVariationalFormulation
;
...
...
@@ -106,7 +107,8 @@ namespace HappyHeart
*
* For instance new velocity -> current velocity.
*/
void
PrepareNextTimeIteration
(
const
ExplicitStepVariationalFormulation
&
varf
);
template
<
class
HyperelasticLawT
>
void
PrepareNextTimeIteration
(
const
ExplicitStepVariationalFormulation
<
HyperelasticLawT
>&
varf
);
public:
...
...
Sources/ModelInstances/UnderDevelopment/Poromechanics/Data/FluidVelocity.hxx
View file @
e4f12a02
...
...
@@ -98,6 +98,24 @@ namespace HappyHeart
{
return
const_cast
<
GlobalVector
&>
(
GetVelFtr
());
}
template
<
class
HyperelasticLawT
>
void
FluidVelocity
::
PrepareNextTimeIteration
(
const
ExplicitStepVariationalFormulation
<
HyperelasticLawT
>&
explicit_varf
)
{
auto
&
current
=
GetNonCstCurrent
();
// \todo #820 Swap...
current
.
Copy
(
GetNew
(),
__FILE__
,
__LINE__
);
Wrappers
::
Petsc
::
MatMult
(
explicit_varf
.
GetMassMatrixWithPorosity
(),
current
,
GetNonCstMassTimesCurrent
(),
__FILE__
,
__LINE__
);
}
...
...
Sources/ModelInstances/UnderDevelopment/Poromechanics/ExplicitStep/VariationalFormulation.cpp
deleted
100644 → 0
View file @
518301ad
//! \file
//
//
// ExplicitStepVariationalFormulation.cpp
// HappyHeart
//
// Created by Sebastien Gilles on 10/03/16.
// Copyright © 2016 Inria. All rights reserved.
//
#include
"FiniteElement/BoundaryConditions/DirichletBoundaryConditionManager.hpp"
#include
"ParameterInstances/GradientBasedElasticityTensor/Internal/Configuration.hpp"
#include
"ModelInstances/UnderDevelopment/Poromechanics/ExplicitStep/VariationalFormulation.hpp"
namespace
HappyHeart
{
namespace
PoromechanicsNS
{
ExplicitStepVariationalFormulation
::
ExplicitStepVariationalFormulation
(
const
Wrappers
::
Mpi
&
mpi
,
const
NumberingSubset
&
numbering_subset
,
const
TimeManager
&
time_manager
,
const
GodOfDof
&
god_of_dof
,
DirichletBoundaryCondition
::
vector_shared_ptr
&&
boundary_condition_list
,
::
HappyHeart
::
PoromechanicsNS
::
VariableHolder
&
variable_holder
,
const
DataNS
::
Fluidmass
&
fluid_mass_data
,
const
DataNS
::
SolidVelocity
&
solid_velocity_data
,
const
DataNS
::
Porosity
&
porosity_data
,
const
DataNS
::
FluidPressure
&
fluid_pressure_data
,
DataNS
::
FluidVelocity
&
fluid_velocity_data
)
:
parent
(
mpi
,
time_manager
,
god_of_dof
,
std
::
move
(
boundary_condition_list
)),
variable_holder_parent
(
variable_holder
),
fluid_mass_data_parent
(
fluid_mass_data
),
fluid_velocity_data_parent
(
fluid_velocity_data
),
solid_velocity_data_parent
(
solid_velocity_data
),
fluid_pressure_data_parent
(
fluid_pressure_data
),
porosity_data_parent
(
porosity_data
),
numbering_subset_
(
numbering_subset
)
{
}
const
GlobalVector
&
ExplicitStepVariationalFormulation
::
Perform
()
{
auto
&
ret
=
GetNonCstExplicitContribution
();
decltype
(
auto
)
variable_holder
=
this
->
GetVariableHolder
();
decltype
(
auto
)
numbering_subset
=
GetNumberingSubset
();
UpdateMatricesAndVectors
();
ApplyEssentialBoundaryCondition
<
VariationalFormulationNS
::
On
::
system_matrix_and_rhs
,
BoundaryConditionMethod
::
penalization
>
(
numbering_subset
,
numbering_subset
);
decltype
(
auto
)
system_matrix
=
parent
::
GetNonCstSystemMatrix
(
numbering_subset
,
numbering_subset
);
system_matrix
.
Assembly
(
__FILE__
,
__LINE__
);
variable_holder
.
DevPrint
<
DevPhase
::
none
>
(
system_matrix
,
"explicit_step_matrix"
);
{
// I now "modify" the system matrix to compute what is called 'ExplicitStepMatrixDirichletNH' in Freefem.
// \todo #820 As noted above, it's likely the same outcome could be written more properly with inohomogeneous
// Dirichlet boundary condition.
decltype
(
auto
)
bc_manager
=
DirichletBoundaryConditionManager
::
GetInstance
();
decltype
(
auto
)
bc
=
bc_manager
.
GetDirichletBoundaryCondition
(
EnumUnderlyingType
(
BoundaryConditionIndex
::
fluid_robin_interface
));
parent
::
GetGodOfDof
().
template
ApplyBoundaryCondition
<
BoundaryConditionMethod
::
penalization
>(
bc
,
system_matrix
);
system_matrix
.
Assembly
(
__FILE__
,
__LINE__
);
variable_holder
.
DevPrint
<
DevPhase
::
none
>
(
system_matrix
,
"explicit_step_matrix_dirichlet_NH"
);
}
SolveLinear
<
IsFactorized
::
no
>
(
numbering_subset
,
numbering_subset
,
__FILE__
,
__LINE__
);
auto
&
solution
=
parent
::
GetNonCstSystemSolution
(
numbering_subset
);
decltype
(
auto
)
fluid_velocity_data
=
GetNonCstFluidVelocityData
();
decltype
(
auto
)
velFtr
=
fluid_velocity_data
.
GetNonCstVelFtr
();
velFtr
.
Copy
(
solution
,
__FILE__
,
__LINE__
);
// \todo #820 Swap!
variable_holder
.
DevPrint
<
DevPhase
::
none
>
(
solution
,
"explicit_step_solution_part_1"
);
// Outcome of the solve is not what we actually need; other contributions must be added to this output.
// Rather than defining a linear operator as in Freefem, we use the already known bilinear one.
decltype
(
auto
)
solid_velocity_data
=
GetSolidVelocityData
();
Wrappers
::
Petsc
::
AXPY
(
1.
,
solid_velocity_data
.
GetHalfSumOnFluidZeroedOutOfRobinInterface
(),
velFtr
,
__FILE__
,
__LINE__
);
variable_holder
.
DevPrint
<
DevPhase
::
none
>
(
velFtr
,
"explicit_step_solution_part_2"
);
Wrappers
::
Petsc
::
MatMult
(
GetSystemMatrixBeforeApplyingBoundaryConditions
(),
velFtr
,
ret
,
__FILE__
,
__LINE__
);
ret
.
Scale
(
-
1.
,
__FILE__
,
__LINE__
);
variable_holder
.
DevPrint
<
DevPhase
::
none
>
(
ret
,
"vector_from_explicit_without_fluidOldVelocity"
);
Wrappers
::
Petsc
::
AXPY
(
1.
,
GetFluidVelocityData
().
GetMassTimesCurrent
(),
ret
,
__FILE__
,
__LINE__
);
variable_holder
.
DevPrint
<
DevPhase
::
none
>
(
ret
,
"vector_from_explicit"
);
return
ret
;
}
void
ExplicitStepVariationalFormulation
::
SupplInit
(
const
InputParameterList
&
input_parameter_data
)
{
decltype
(
auto
)
god_of_dof
=
parent
::
GetGodOfDof
();
decltype
(
auto
)
mesh
=
god_of_dof
.
GetGeometricMeshRegion
();
decltype
(
auto
)
variable_holder
=
GetVariableHolder
();
{
using
parameter_type
=
::
HappyHeart
::
Internal
::
ParameterNS
::
ParameterInstance
<
ParameterNS
::
Type
::
scalar
,
ParameterNS
::
Policy
::
Constant
,
ParameterNS
::
TimeDependencyNS
::
None
>
;
pseudo_young_modulus_
=
std
::
make_unique
<
parameter_type
>
(
"Pseudo young modulus"
,
mesh
,
2.
*
variable_holder
.
GetFluidViscosity
().
GetConstantValue
());
pseudo_poisson_ratio_
=
std
::
make_unique
<
parameter_type
>
(
"Zero poisson ratio"
,
mesh
,
0.
);
}
namespace
ipl
=
Utilities
::
InputParameterListNS
;
fluid_source_term_
=
ipl
::
Extract
<
InputParameter
::
PoromechanicsNS
::
FluidSourceTerm
>::
Value
(
input_parameter_data
);
DefineOperators
();
}
void
ExplicitStepVariationalFormulation
::
UpdateMatricesAndVectors
()
{
decltype
(
auto
)
numbering_subset
=
GetNumberingSubset
();
auto
&
system_matrix
=
parent
::
GetNonCstSystemMatrix
(
numbering_subset
,
numbering_subset
);
auto
&
system_rhs
=
parent
::
GetNonCstSystemRhs
(
numbering_subset
);
decltype
(
auto
)
time_manager
=
parent
::
GetTimeManager
();
decltype
(
auto
)
variable_holder
=
this
->
GetVariableHolder
();
const
auto
time_step
=
time_manager
.
GetTimeStep
();
system_matrix
.
ZeroEntries
(
__FILE__
,
__LINE__
);
system_rhs
.
ZeroEntries
(
__FILE__
,
__LINE__
);
{
auto
&
mass_matrix
=
GetNonCstMassMatrixWithPorosity
();
mass_matrix
.
ZeroEntries
(
__FILE__
,
__LINE__
);
assert
(
!
NumericNS
::
IsZero
(
time_step
));
const
double
coefficient
=
variable_holder
.
GetFluidDensity
().
GetConstantValue
()
/
time_step
;
GlobalMatrixWithCoefficient
with_coeff_matrix
(
mass_matrix
,
coefficient
);
GetMassOperatorWithPorosity
().
Assemble
(
std
::
make_tuple
(
std
::
ref
(
with_coeff_matrix
)));
variable_holder
.
DevPrint
<
DevPhase
::
none
>
(
mass_matrix
,
"explicit_step_mass_matrix"
);
system_matrix
.
Copy
(
mass_matrix
,
__FILE__
,
__LINE__
);
}
{
auto
&
mass_matrix
=
GetNonCstMassMatrixWithPressure
();
mass_matrix
.
ZeroEntries
(
__FILE__
,
__LINE__
);
const
double
coefficient
=
variable_holder
.
GetFluidDensity
().
GetConstantValue
()
*
GetFluidSourceTerm
();
GlobalMatrixWithCoefficient
with_coeff_matrix
(
mass_matrix
,
coefficient
);
GetMassOperatorWithPressure
().
Assemble
(
std
::
make_tuple
(
std
::
ref
(
with_coeff_matrix
)));
// variable_holder.DevPrint<DevPhase::none>(mass_matrix,
// "explicit_step_mass_matrix");
Wrappers
::
Petsc
::
AXPY
<
NonZeroPattern
::
same
>
(
1.
,
mass_matrix
,
system_matrix
,
__FILE__
,
__LINE__
);
}
{
GlobalMatrixWithCoefficient
with_coeff_matrix
(
system_matrix
,
1.
);
GetViscosityOperator
().
Assemble
(
std
::
make_tuple
(
std
::
ref
(
with_coeff_matrix
)));
}
decltype
(
auto
)
fluid_velocity_data
=
this
->
GetFluidVelocityData
();
decltype
(
auto
)
solid_velocity_data
=
this
->
GetSolidVelocityData
();
// decltype(auto) fluid_pressure_data = this->GetFluidPressureData();
variable_holder
.
DevPrint
<
DevPhase
::
none
>
(
system_matrix
,
"explicit_step_matrix_no_ale_without_beta_contrib"
);
variable_holder
.
DevPrint
<
DevPhase
::
none
>
(
system_matrix
,
"explicit_step_matrix_no_ale_with_beta_contrib"
);
{
GlobalMatrixWithCoefficient
with_coeff_matrix
(
system_matrix
,
1.
);
decltype
(
auto
)
current_velocity
=
fluid_velocity_data
.
GetCurrent
();
variable_holder
.
DevPrint
<
DevPhase
::
none
>
(
current_velocity
,
"explicit_step_prev_time_it_fluid_velocity"
);
variable_holder
.
DevPrint
<
DevPhase
::
none
>
(
solid_velocity_data
.
GetHalfSumOnFluid
(),
"explicit_step_velSHalfVhfr"
);
// variable_holder.DevPrint<DevPhase::none>(fluid_pressure_data.GetCurrent(),
// "explicit_step_pressure_prev_time_it");
auto
&
modified_velocity
=
GetNonCstModifiedCurrentFluidVelocity
();
modified_velocity
.
Copy
(
current_velocity
,
__FILE__
,
__LINE__
);
Wrappers
::
Petsc
::
AXPY
(
-
1.
,
fluid_velocity_data
.
GetHalfSum
(),
modified_velocity
,
__FILE__
,
__LINE__
);
GetAleOperator
().
Assemble
(
std
::
make_tuple
(
std
::
ref
(
with_coeff_matrix
)),
modified_velocity
);
}
{
// In Freefem: explicitStepVector += fluidOldVelocity;
Wrappers
::
Petsc
::
AXPY
(
1.
,
fluid_velocity_data
.
GetMassTimesCurrent
(),
system_rhs
,
__FILE__
,
__LINE__
);
variable_holder
.
DevPrint
<
DevPhase
::
none
>
(
system_rhs
,
"explicit_step_vector_part_2"
);
}
{
// \todo #820 It's highly possible that it can be done more elegantly with a variable Dirichlet condition
// (it's possible Bruno did current scheme for conveniency with Freefem interface).
// \todo I do not update the value of GetHalfSumOnFluidZeroedOutOfRobinInterface() here; if it's necessary
// just call:
// solid_velocity_data.UpdateHalfSum();
decltype
(
auto
)
prescribed_vel
=
solid_velocity_data
.
GetHalfSumOnFluidZeroedOutOfRobinInterface
();
variable_holder
.
DevPrint
<
DevPhase
::
none
>
(
prescribed_vel
,
"explicit_step_prescribed_velocity"
);
DataNS
::
FluidVelocity
::
global_vector_temporary_manager
helper
(
fluid_velocity_data
);
auto
&
modified_prescribed_vel
=
helper
.
GetNonCstVector
();
Wrappers
::
Petsc
::
MatMult
(
system_matrix
,
prescribed_vel
,
modified_prescribed_vel
,
__FILE__
,
__LINE__
);
variable_holder
.
DevPrint
<
DevPhase
::
none
>
(
modified_prescribed_vel
,
"explicit_step_modified_prescribed_velocity"
);
Wrappers
::
Petsc
::
AXPY
(
-
1.
,
modified_prescribed_vel
,
system_rhs
,
__FILE__
,
__LINE__
);
}
variable_holder
.
DevPrint
<
DevPhase
::
none
>
(
system_rhs
,
"explicit_step_rhs"
);
GetNonCstSystemMatrixBeforeApplyingBoundaryConditions
().
Copy
(
system_matrix
,
__FILE__
,
__LINE__
);
}
void
ExplicitStepVariationalFormulation
::
DefineOperators
()
{
decltype
(
auto
)
god_of_dof
=
parent
::
GetGodOfDof
();
const
auto
&
felt_space
=
god_of_dof
.
GetFEltSpace
(
EnumUnderlyingType
(
FEltSpaceIndex
::
fluid
));
decltype
(
auto
)
unknown_manager
=
UnknownManager
::
GetInstance
();
decltype
(
auto
)
fluid_velocity
=
unknown_manager
.
GetUnknown
(
EnumUnderlyingType
(
UnknownIndex
::
fluid_velocity
));
decltype
(
auto
)
mesh
=
god_of_dof
.
GetGeometricMeshRegion
();
const
auto
mesh_dimension
=
mesh
.
GetDimension
();
decltype
(
auto
)
porosity_data
=
this
->
GetPorosityData
();
decltype
(
auto
)
porosity
=
porosity_data
.
GetNewOnFluidAsParam
();
decltype
(
auto
)
variable_holder
=
GetVariableHolder
();
decltype
(
auto
)
fluid_density
=
variable_holder
.
GetFluidDensity
();
decltype
(
auto
)
fluid_pressure_data
=
this
->
GetFluidPressureData
();
mass_operator_with_porosity_
=
std
::
make_unique
<
mass_with_porosity_operator_type
>
(
felt_space
,
fluid_velocity
,
porosity
);
mass_operator_with_pressure_
=
std
::
make_unique
<
mass_with_pressure_operator_type
>
(
felt_space
,
fluid_velocity
,
fluid_pressure_data
.
GetCurrentAsParam
());
using
op_type
=
GlobalVariationalOperatorNS
::
GradOnGradientBasedElasticityTensor
;
const
auto
configuration
=
mesh_dimension
==
2
?
ParameterNS
::
GradientBasedElasticityTensorConfiguration
::
dim2_plane_strain
:
ParameterNS
::
GradientBasedElasticityTensorConfiguration
::
dim3
;
viscosity_operator_
=
std
::
make_unique
<
op_type
>
(
felt_space
,
fluid_velocity
,
porosity
,
GetPseudoYoungModulus
(),
GetPseudoPoissonRatio
(),
configuration
);
ale_operator_
=
std
::
make_unique
<
GlobalVariationalOperatorNS
::
Ale
>
(
felt_space
,
fluid_velocity
,
porosity_data
.
GetCurrentOnFluidAsParam
(),
fluid_density
);
}
void
ExplicitStepVariationalFormulation
::
AllocateMatricesAndVectors
()
{
const
auto
&
numbering_subset
=
GetNumberingSubset
();
parent
::
AllocateSystemMatrix
(
numbering_subset
,
numbering_subset
);
parent
::
AllocateSystemVector
(
numbering_subset
);
decltype
(
auto
)
system_matrix
=
parent
::
GetSystemMatrix
(
numbering_subset
,
numbering_subset
);
mass_matrix_with_porosity_
=
std
::
make_unique
<
GlobalMatrix
>
(
system_matrix
);
system_matrix_before_applying_bc_
=
std
::
make_unique
<
GlobalMatrix
>
(
system_matrix
);
mass_matrix_with_pressure_
=
std
::
make_unique
<
GlobalMatrix
>
(
system_matrix
);
decltype
(
auto
)
system_rhs
=
parent
::
GetSystemRhs
(
numbering_subset
);
modified_fluid_current_velocity_
=
std
::
make_unique
<
GlobalVector
>
(
system_rhs
);
explicit_contribution_
=
std
::
make_unique
<
GlobalVector
>
(
system_rhs
);
}
}
// namespace PoromechanicsNS
}
// namespace HappyHeart
Sources/ModelInstances/UnderDevelopment/Poromechanics/ExplicitStep/VariationalFormulation.hpp
View file @
e4f12a02
...
...
@@ -18,7 +18,10 @@
# include "Geometry/Domain/Domain.hpp"
# include "FiniteElement/BoundaryConditions/DirichletBoundaryConditionManager.hpp"
# include "OperatorInstances/VariationalOperator/BilinearForm/VariableMass.hpp"
# include "ParameterInstances/GradientBasedElasticityTensor/Internal/Configuration.hpp"
# include "FormulationSolver/VariationalFormulation.hpp"
...
...
@@ -32,6 +35,7 @@
# include "ModelInstances/UnderDevelopment/Poromechanics/Crtp/SolidVelocityData.hpp"
# include "ModelInstances/UnderDevelopment/Poromechanics/Crtp/PorosityData.hpp"
# include "ModelInstances/UnderDevelopment/Poromechanics/Crtp/FluidPressureData.hpp"
# include "ModelInstances/UnderDevelopment/Poromechanics/Crtp/NewFluidPressureData.hpp"
namespace
HappyHeart
...
...
@@ -42,43 +46,53 @@ namespace HappyHeart
{
template
<
class
HyperelasticLawT
>
class
ExplicitStepVariationalFormulation
:
public
HappyHeart
::
VariationalFormulation
<
ExplicitStepVariationalFormulation
,
ExplicitStepVariationalFormulation
<
HyperelasticLawT
>
,
EnumUnderlyingType
(
SolverIndex
::
direct
)
>
,
public
::
HappyHeart
::
PoromechanicsNS
::
Crtp
::
VariableHolder
<
ExplicitStepVariationalFormulation
>
,
public
::
HappyHeart
::
PoromechanicsNS
::
Crtp
::
VariableHolder
<
ExplicitStepVariationalFormulation
<
HyperelasticLawT
>
>
,
public
::
HappyHeart
::
PoromechanicsNS
::
Crtp
::
FluidmassData
<
ExplicitStepVariationalFormulation
,
ExplicitStepVariationalFormulation
<
HyperelasticLawT
>
,
Utilities
::
Access
::
read_only
>
,
public
::
HappyHeart
::
PoromechanicsNS
::
Crtp
::
FluidVelocityData
<
ExplicitStepVariationalFormulation
,
ExplicitStepVariationalFormulation
<
HyperelasticLawT
>
,
Utilities
::
Access
::
read_and_write
>
,
public
::
HappyHeart
::
PoromechanicsNS
::
Crtp
::
SolidVelocityData
<
ExplicitStepVariationalFormulation
,
ExplicitStepVariationalFormulation
<
HyperelasticLawT
>
,
Utilities
::
Access
::
read_only
>
,
public
::
HappyHeart
::
PoromechanicsNS
::
Crtp
::
FluidPressureData
<
ExplicitStepVariationalFormulation
,
ExplicitStepVariationalFormulation
<
HyperelasticLawT
>
,
Utilities
::
Access
::
read_only
>
,