Commit b9638375 authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#1049 Fix Fiber and bidomain model.

parent c361a3db
......@@ -156,6 +156,40 @@ Domain1 = {
} -- Domain1
Domain2 = {
-- Index of the geometric mesh upon which the domain is defined (as defined in the present file). Might be
-- left empty if domain not limited to one mesh; at most one value is expected here.
-- Expected format: {VALUE1, VALUE2, ...}
mesh_index = { 1 },
-- List of dimensions encompassed by the domain. Might be left empty if no restriction at all upon
-- dimensions.
-- Expected format: {VALUE1, VALUE2, ...}
-- Constraint: ops_in(v, {0, 1, 2, 3})
dimension_list = { },
-- List of mesh labels encompassed by the domain. Might be left empty if no restriction at all upon mesh
-- labels. This parameter does not make sense if no mesh is defined for the domain.
-- Expected format: {VALUE1, VALUE2, ...}
mesh_label_list = { },
-- List of geometric element types considered in the domain. Might be left empty if no restriction upon
-- these. No constraint is applied at Ops level, as some geometric element types could be added after
-- generation of current input parameter file. Current list is below; if an incorrect value is put there it
-- will be detected a bit later when the domain object is built.
-- The known types when this file was generated are:
-- . Point1
-- . Segment2, Segment3
-- . Triangle3, Triangle6
-- . Quadrangle4, Quadrangle8, Quadrangle9
-- . Tetrahedron4, Tetrahedron10
-- . Hexahedron8, Hexahedron20, Hexahedron27.
-- Expected format: {"VALUE1", "VALUE2", ...}
geometric_element_type_list = { }
} -- Domain2
FiniteElementSpace1 = {
-- Index of the god of dof into which the finite element space is defined.
......@@ -436,9 +470,9 @@ Fiber_vector_1 = {
-- Expected format: "VALUE"
ensight_file = "${HOME}/Codes/HappyHeart/Data/Fiber/square_tria_4kv_fibers_at_zero.vct",
-- Index of the geometric mesh region upon which parameter is defined.
-- Index of the domain upon which parameter is defined.
-- Expected format: VALUE
mesh_index = 1,
domain_index = 2,
-- Index of the finite element space upon which parameter is defined.
-- Expected format: VALUE
......
......@@ -156,6 +156,42 @@ Domain1 = {
} -- Domain1
Domain2 = {
-- Index of the geometric mesh upon which the domain is defined (as defined in the present file). Might be
-- left empty if domain not limited to one mesh; at most one value is expected here.
-- Expected format: {VALUE1, VALUE2, ...}
mesh_index = { 1 },
-- List of dimensions encompassed by the domain. Might be left empty if no restriction at all upon
-- dimensions.
-- Expected format: {VALUE1, VALUE2, ...}
-- Constraint: ops_in(v, {0, 1, 2, 3})
dimension_list = { },
-- List of mesh labels encompassed by the domain. Might be left empty if no restriction at all upon mesh
-- labels. This parameter does not make sense if no mesh is defined for the domain.
-- Expected format: {VALUE1, VALUE2, ...}
mesh_label_list = { },
-- List of geometric element types considered in the domain. Might be left empty if no restriction upon
-- these. No constraint is applied at Ops level, as some geometric element types could be added after
-- generation of current input parameter file. Current list is below; if an incorrect value is put there it
-- will be detected a bit later when the domain object is built.
-- The known types when this file was generated are:
-- . Point1
-- . Segment2, Segment3
-- . Triangle3, Triangle6
-- . Quadrangle4, Quadrangle8, Quadrangle9
-- . Tetrahedron4, Tetrahedron10
-- . Hexahedron8, Hexahedron20, Hexahedron27.
-- Expected format: {"VALUE1", "VALUE2", ...}
geometric_element_type_list = { }
} -- Domain2
FiniteElementSpace1 = {
-- Index of the god of dof into which the finite element space is defined.
......@@ -434,7 +470,7 @@ Fiber_vector_1 = {
-- Expected format: "VALUE"
ensight_file = "${HOME}/Codes/HappyHeart/Data/Fiber/square_tria_4kv_fibers_pi_4.vct",
-- Index of the geometric mesh region upon which parameter is defined.
-- Index of the domain upon which parameter is defined.
-- Expected format: VALUE
mesh_index = 1,
......
......@@ -730,7 +730,7 @@ Fiber_vector_1 = {
-- Expected format: "VALUE"
ensight_file = "${HOME}/Codes/HappyHeart/Data/Fiber/Zygot_Vent_TeleSystol_6Tet_SurfElecSimu_VolVRVL_cmFibers_refbeforeInterp.00000.vct",
-- Index of the geometric mesh region upon which parameter is defined.
-- Index of the domain upon which parameter is defined.
-- Expected format: VALUE
mesh_index = 1,
......
......@@ -759,7 +759,7 @@ Fiber_vector_1 = {
-- Expected format: "VALUE"
ensight_file = "${HOME}/Codes/HappyHeart/Data/Fiber/ventriVol300Fib.vct",
-- Index of the geometric mesh region upon which parameter is defined.
-- Index of the domain upon which parameter is defined.
-- Expected format: VALUE
mesh_index = 1,
......
......@@ -707,7 +707,7 @@ Fiber_vector_1 = {
-- Expected format: "VALUE"
ensight_file = "${HOME}/Codes/HappyHeart/Data/Fiber/Zygot_AtriaDbl_TeleSystol_8TetInVent_SurfSimuElec_fibers.00000.vct",
-- Index of the geometric mesh region upon which parameter is defined.
-- Index of the domain upon which parameter is defined.
-- Expected format: VALUE
mesh_index = 1,
......@@ -728,7 +728,7 @@ Fiber_scalar_2 = {
-- Expected format: "VALUE"
ensight_file = "${HOME}/Codes/HappyHeart/Data/Fiber/Zygot_AtriaDbl_TeleSystol_8TetInVent_SurfSimuElec_angles.00000.scl",
-- Index of the geometric mesh region upon which parameter is defined.
-- Index of the domain upon which parameter is defined.
-- Expected format: VALUE
mesh_index = 1,
......
......@@ -707,7 +707,7 @@ Fiber_vector_1 = {
-- Expected format: "VALUE"
ensight_file = "${HOME}/Codes/HappyHeart/Data/Fiber/Zygot_AtriaDbl_TeleSystol_8TetInVent_SurfSimuElec_fibers.00000.vct",
-- Index of the geometric mesh region upon which parameter is defined.
-- Index of the domain upon which parameter is defined.
-- Expected format: VALUE
mesh_index = 1,
......@@ -728,7 +728,7 @@ Fiber_scalar_2 = {
-- Expected format: "VALUE"
ensight_file = "${HOME}/Codes/HappyHeart/Data/Fiber/Zygot_AtriaDbl_TeleSystol_8TetInVent_SurfSimuElec_angles.00000.scl",
-- Index of the geometric mesh region upon which parameter is defined.
-- Index of the domain upon which parameter is defined.
-- Expected format: VALUE
mesh_index = 1,
......
......@@ -170,7 +170,7 @@ Fiber_vector_1 = {
ensight_file = '${HOME}/Codes/HappyHeart/Data/Fiber/fibers.vct',
--ensight_file = "/Volumes/Data/sebastien/Fiber/Zygot_AtriaDbl_TeleSystol_8TetInVent_SurfSimuElec_fibers.00000.vct",
-- Index of the geometric mesh region upon which parameter is defined.
-- Index of the domain upon which parameter is defined.
-- Expected format: VALUE
mesh_index = 1,
......@@ -192,7 +192,7 @@ Fiber_scalar_2 = {
ensight_file = '${HOME}/Codes/HappyHeart/Data/Fiber/distance.scl',
--ensight_file = "/Volumes/Data/sebastien/Fiber/Zygot_AtriaDbl_TeleSystol_8TetInVent_SurfSimuElec_distance.00000.scl",
-- Index of the geometric mesh region upon which parameter is defined.
-- Index of the domain upon which parameter is defined.
-- Expected format: VALUE
mesh_index = 1,
......
......@@ -82,12 +82,12 @@ namespace HappyHeart
/*!
* \brief Index of the geometric mesh region onto which fiber is defined. It is expected this parameter
* \brief Index of the \a Domain onto which fiber is defined. It is expected this parameter
* is compatible with the fiber file.
*/
struct GeometricMeshRegionIndex
: public Utilities::InputParameterListNS::Crtp::InputParameter<GeometricMeshRegionIndex, self, FiberNS::GeometricMeshRegionIndex::storage_type>,
public FiberNS::GeometricMeshRegionIndex
struct DomainIndex
: public Utilities::InputParameterListNS::Crtp::InputParameter<DomainIndex, self, FiberNS::DomainIndex::storage_type>,
public FiberNS::DomainIndex
{ };
......@@ -117,7 +117,7 @@ namespace HappyHeart
<
EnsightFile,
// Nature,
GeometricMeshRegionIndex,
DomainIndex,
FEltSpaceIndex,
UnknownName
>;
......
......@@ -53,26 +53,26 @@ namespace HappyHeart
}
const std::string& GeometricMeshRegionIndex::NameInFile()
const std::string& DomainIndex::NameInFile()
{
static std::string ret("mesh_index");
static std::string ret("domain_index");
return ret;
}
const std::string& GeometricMeshRegionIndex::Description()
const std::string& DomainIndex::Description()
{
static std::string ret("Index of the geometric mesh region upon which parameter is defined." );
static std::string ret("Index of the domain upon which parameter is defined." );
return ret;
}
const std::string& GeometricMeshRegionIndex::Constraint()
const std::string& DomainIndex::Constraint()
{
return Utilities::EmptyString();
}
const std::string& GeometricMeshRegionIndex::DefaultValue()
const std::string& DomainIndex::DefaultValue()
{
return Utilities::EmptyString();
}
......
......@@ -80,7 +80,7 @@ namespace HappyHeart
* \brief Index of the geometric mesh region.
*
*/
struct GeometricMeshRegionIndex
struct DomainIndex
{
using storage_type = unsigned int;
......
......@@ -35,7 +35,6 @@ namespace HappyHeart
{ }
template <Advanced::ReactionLawNS::ReactionLawName ReactionLawNameT>
void BidomainVariationalFormulation<ReactionLawNameT>::AllocateMatricesAndVectors()
{
......@@ -60,32 +59,30 @@ namespace HappyHeart
template <Advanced::ReactionLawNS::ReactionLawName ReactionLawNameT>
void BidomainVariationalFormulation<ReactionLawNameT>::SupplInit(const InputParameterList& input_parameter_data)
{
const auto& god_of_dof = this->GetGodOfDof();
const auto& mesh = god_of_dof.GetGeometricMeshRegion();
decltype(auto) domain = DomainManager::GetInstance().GetDomain(EnumUnderlyingType(DomainIndex::full_mesh));
using Diffusion = InputParameter::Diffusion;
intracellular_trans_diffusion_tensor_ =
InitScalarParameterFromInputData<Diffusion::Tensor<EnumUnderlyingType(TensorIndex::intracellular_trans_diffusion_tensor)>>("Intracellular Trans Diffusion tensor",
mesh,
input_parameter_data);
domain,
input_parameter_data);
extracellular_trans_diffusion_tensor_ =
InitScalarParameterFromInputData<Diffusion::Tensor<EnumUnderlyingType(TensorIndex::extracellular_trans_diffusion_tensor)>>("Extracellular Trans Diffusion tensor",
mesh,
input_parameter_data);
domain,
input_parameter_data);
intracellular_fiber_diffusion_tensor_ =
InitScalarParameterFromInputData<Diffusion::Tensor<EnumUnderlyingType(TensorIndex::intracellular_fiber_diffusion_tensor)>>("Intracellular Fiber Diffusion tensor",
mesh,
input_parameter_data);
domain,
input_parameter_data);
extracellular_fiber_diffusion_tensor_ =
InitScalarParameterFromInputData<Diffusion::Tensor<EnumUnderlyingType(TensorIndex::extracellular_fiber_diffusion_tensor)>>("Extracellular Fiber Diffusion tensor",
mesh,
input_parameter_data);
domain,
input_parameter_data);
transcellular_density_ =
InitScalarParameterFromInputData<Diffusion::Density>("Transcellular Density",
mesh,
input_parameter_data);
domain,
input_parameter_data);
if (!GetTranscellularDiffusionDensity().IsConstant())
throw Exception("Current Bidomain model is restricted to a constant diffusion density.",
......@@ -131,7 +128,9 @@ namespace HappyHeart
void BidomainVariationalFormulation<ReactionLawNameT>::DefineOperators(const InputParameterList& input_parameter_data)
{
const auto& god_of_dof = this->GetGodOfDof();
const auto& mesh = god_of_dof.GetGeometricMeshRegion();
decltype(auto) domain = DomainManager::GetInstance().GetDomain(EnumUnderlyingType(DomainIndex::full_mesh));
const auto& felt_space_highest_dimension = god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::highest_dimension));
const auto& transcellular_potential = UnknownManager::GetInstance().GetUnknown(EnumUnderlyingType(UnknownIndex::transcellular_potential));
......@@ -163,7 +162,7 @@ namespace HappyHeart
extracellular_potential);
reaction_law_ = std::make_unique<reaction_law_type>(input_parameter_data,
mesh,
domain,
this->GetTimeManager(),
felt_space_highest_dimension.GetQuadratureRulePerTopology());
......@@ -175,8 +174,8 @@ namespace HappyHeart
using parameter_type_on_square = InputParameter::TransientSource<EnumUnderlyingType(ForceIndexList::transcellular_current_applied_on_square)>;
transcellular_current_applied_on_square_parameter_ =
InitThreeDimensionalParameterFromInputData<parameter_type_on_square>("Current Applied On Square",
mesh,
input_parameter_data);
domain,
input_parameter_data);
if (transcellular_current_applied_on_square_parameter_ != nullptr)
{
......
......@@ -50,7 +50,8 @@ namespace HappyHeart
enum class DomainIndex
{
highest_dimension = 1
highest_dimension = 1,
full_mesh = 2
};
......@@ -124,6 +125,7 @@ namespace HappyHeart
InputParameter::Mesh<EnumUnderlyingType(MeshIndex::mesh)>,
InputParameter::Domain<EnumUnderlyingType(DomainIndex::highest_dimension)>,
InputParameter::Domain<EnumUnderlyingType(DomainIndex::full_mesh)>,
InputParameter::FEltSpace<EnumUnderlyingType(FEltSpaceIndex::highest_dimension)>,
......
......@@ -143,7 +143,7 @@ namespace HappyHeart
*/
explicit FiberList(unsigned int unique_id,
const std::string& fiber_file,
const GeometricMeshRegion& geometric_mesh_region,
const Domain& domain,
const FEltSpace& felt_space,
const TimeManager& time_manager,
const Unknown& unknown);
......
......@@ -27,13 +27,13 @@ namespace HappyHeart
template<ParameterNS::Type TypeT>
FiberList<TypeT>::FiberList(const unsigned int unique_id,
const std::string& fiber_file,
const GeometricMeshRegion& geometric_mesh_region,
const Domain& domain,
const FEltSpace& felt_space,
const TimeManager& time_manager,
const Unknown& unknown)
: unique_id_parent(unique_id),
parent("Fiber manager",
geometric_mesh_region)
domain)
{
const auto god_of_dof_ptr = felt_space.GetGodOfDofFromWeakPtr();
const auto& god_of_dof = *god_of_dof_ptr;
......@@ -43,9 +43,11 @@ namespace HappyHeart
global_vector_ = std::make_unique<GlobalVector>(numbering_subset);
AllocateGlobalVector(god_of_dof, *global_vector_);
assert(god_of_dof.GetGeometricMeshRegion().GetUniqueId() == geometric_mesh_region.GetUniqueId());
assert(god_of_dof.GetGeometricMeshRegion().GetUniqueId() == domain.GetGeometricMeshRegion().GetUniqueId());
typename Internal::FiberNS::Traits<TypeT>::value_list_per_coord_index_type value_list_per_coord_index;
decltype(auto) geometric_mesh_region = god_of_dof.GetGeometricMeshRegion();
Internal::FiberNS::ReadFiberFile<TypeT>(god_of_dof.MpiHappyHeart(),
fiber_file,
......@@ -64,7 +66,7 @@ namespace HappyHeart
global_vector);
ParameterAtDof<TypeT, ParameterNS::TimeDependencyNS::None> param_at_dof("Temporary",
geometric_mesh_region,
domain,
felt_space,
unknown,
global_vector);
......@@ -78,7 +80,7 @@ namespace HappyHeart
>;
underlying_parameter_ = converter_type::Perform("Fiber manager",
geometric_mesh_region,
domain,
felt_space.GetQuadratureRulePerTopology(),
time_manager,
param_at_dof);
......
......@@ -95,7 +95,7 @@ namespace HappyHeart
*/
void Create(const unsigned int unique_id,
const std::string& fiber_file,
const GeometricMeshRegion& geometric_mesh_region,
const Domain& domain,
const FEltSpace& felt_space,
const Unknown& unknown);
......
......@@ -47,19 +47,22 @@ namespace HappyHeart
namespace ipl = Utilities::InputParameterListNS;
decltype(auto) ensight_file = ipl::ExtractPathParameter<typename FiberSectionT::EnsightFile>(section);
decltype(auto) mesh_index = ipl::ExtractParameter<typename FiberSectionT::GeometricMeshRegionIndex>(section);
decltype(auto) domain_index = ipl::ExtractParameter<typename FiberSectionT::DomainIndex>(section);
decltype(auto) felt_space_index = ipl::ExtractParameter<typename FiberSectionT::FEltSpaceIndex>(section);
decltype(auto) unknown_name = ipl::ExtractParameter<typename FiberSectionT::UnknownName>(section);
const auto& mesh = Internal::MeshNS::GeometricMeshRegionManager::GetInstance().GetMesh(mesh_index);
const auto& god_of_dof = GodOfDofManager::GetInstance().GetGodOfDof(mesh_index);
const auto& domain = DomainManager::GetInstance().GetDomain(domain_index);
decltype(auto) mesh = domain.GetGeometricMeshRegion();
const auto& god_of_dof = GodOfDofManager::GetInstance().GetGodOfDof(mesh.GetUniqueId());
const auto& felt_space = god_of_dof.GetFEltSpace(felt_space_index);
const auto& unknown = UnknownManager::GetInstance().GetUnknown(unknown_name);
Create(section.GetUniqueId(),
ensight_file,
mesh,
domain,
felt_space,
unknown);
}
......@@ -68,7 +71,7 @@ namespace HappyHeart
template<ParameterNS::Type TypeT>
void FiberListManager<TypeT>::Create(const unsigned int unique_id,
const std::string& fiber_file,
const GeometricMeshRegion& geometric_mesh_region,
const Domain& domain,
const FEltSpace& felt_space,
const Unknown& unknown)
{
......@@ -76,7 +79,7 @@ namespace HappyHeart
// with both clang and gcc.
fiber_list_type* buf = new fiber_list_type(unique_id,
fiber_file,
geometric_mesh_region,
domain,
felt_space,
GetTimeManager(),
unknown);
......
......@@ -64,8 +64,7 @@ namespace HappyHeart
* tick is used here so it might be a pointer, a reference or a plain chain of characters.
*
* \param[in] name Name given to the parameter, just for output logs.
* \param[in] mesh Mesh upon which both \a Parameter are defined. It is for conveniency: infos is already
* available in \a param_at_dof but not that easily.
* \param[in] domain Domain upon which both \a Parameter are defined.
* \param[in] quadrature_rule_per_topology \a QuadratureRule to use for each topology. Make sure you cover
* all the cases met in the \a FEltSpace underlying \a param_at_dof.
* \copydetails doxygen_hide_time_manager_arg
......@@ -76,7 +75,7 @@ namespace HappyHeart
template<class StringT>
static typename ParameterAtQuadraturePoint<TypeT, TimeDependencyT>::unique_ptr
Perform(StringT&& name,
const GeometricMeshRegion& mesh,
const Domain& domain,
const QuadratureRulePerTopology& quadrature_rule_per_topology,
const TimeManager& time_manager,
const ParameterAtDof<TypeT, TimeDependencyT, NfeltSpaceT>& param_at_dof);
......
......@@ -30,16 +30,18 @@ namespace HappyHeart
typename ParameterAtQuadraturePoint<TypeT, TimeDependencyT>::unique_ptr
FromParameterAtDof<TypeT, TimeDependencyT, NfeltSpaceT>
::Perform(StringT&& name,
const GeometricMeshRegion& mesh,
const Domain& domain,
const QuadratureRulePerTopology& quadrature_rule_per_topology,
const TimeManager& time_manager,
const ParameterAtDof<TypeT, TimeDependencyT, NfeltSpaceT>& param_at_dof)
{
decltype(auto) felt_space_storage = param_at_dof.GetFEltSpaceStorage(); // current class is a friend of \a AtDof policy.
decltype(auto) mesh = domain.GetGeometricMeshRegion();
const auto mesh_dimension = mesh.GetDimension();
auto ret = std::make_unique<ParameterAtQuadraturePoint<TypeT, TimeDependencyT>>(name,
mesh,
domain,
quadrature_rule_per_topology,
Traits<TypeT>::AllocateDefaultValue(mesh_dimension, mesh_dimension),
time_manager);
......
......@@ -15,6 +15,17 @@ MODE='debug'
# Choose either 'static' or 'shared'. 'static' doesn't work for Linux at the moment.
LIBRARY_TYPE='shared'
# Whether specific macros should be activated or not.
# If true, add a (costly) method that gives an hint whether an UpdateGhost() call was relevant or not.
HAPPY_HEART_CHECK_UPDATE_GHOSTS_CALL_RELEVANCE = False
# If true, TimeKeep gains the ability to track times between each call of PrintTimeElapsed(). If not, PrintTimeElapsed() is flatly ignored.
HAPPY_HEART_EXTENDED_TIME_KEEP = False
# If true, there are additional checks that no nan and inf appears in the code. Even if False, solver always check for the validity of its solution (if a nan or an inf is present the SolveLinear() or SolveNonLinear() operation throws with a dedicated Petsc error). Advised in debug mode and up to you in release mode.
HAPPY_HEART_CHECK_NAN_AND_INF = False
# OpenMPI libary.
OPEN_MPI_INCL_DIR='/opt/LibraryVersions/gcc/5.3/Openmpi/include'
OPEN_MPI_LIB_DIR='/opt/LibraryVersions/gcc/5.3/Openmpi/lib'
......
......@@ -21,7 +21,10 @@ LIBRARY_TYPE='static'
HAPPY_HEART_CHECK_UPDATE_GHOSTS_CALL_RELEVANCE = False
# If true, TimeKeep gains the ability to track times between each call of PrintTimeElapsed(). If not, PrintTimeElapsed() is flatly ignored.
HAPPY_HEART_EXTENDED_TIME_KEEP = True
HAPPY_HEART_EXTENDED_TIME_KEEP = False
# If true, there are additional checks that no nan and inf appears in the code. Even if False, solver always check for the validity of its solution (if a nan or an inf is present the SolveLinear() or SolveNonLinear() operation throws with a dedicated Petsc error). Advised in debug mode and up to you in release mode.
HAPPY_HEART_CHECK_NAN_AND_INF = False
# OpenMPI libary.
OPEN_MPI_INCL_DIR='/Users/Shared/LibraryVersions/clang/Openmpi/include'
......
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