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

#1049 Fix ReactionDiffusion model.

parent 54e14c63
......@@ -179,6 +179,39 @@ Domain3 = {
geometric_element_type_list = { }
}
Domain4 = {
-- 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 exxpected 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 = { }
}
-- FiniteElementSpace1
FiniteElementSpace1 = {
-- Index of the god of dof into which the finite element space is defined.
......
......@@ -179,6 +179,40 @@ Domain3 = {
geometric_element_type_list = { }
}
Domain4 = {
-- 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 exxpected 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 = { }
}
-- FiniteElementSpace1
FiniteElementSpace1 = {
-- Index of the god of dof into which the finite element space is defined.
......
......@@ -177,6 +177,40 @@ Domain3 = {
geometric_element_type_list = { }
}
Domain4 = {
-- 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 exxpected 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 = { }
}
-- FiniteElementSpace1
FiniteElementSpace1 = {
-- Index of the god of dof into which the finite element space is defined.
......
......@@ -179,6 +179,39 @@ Domain3 = {
geometric_element_type_list = { }
}
Domain4 = {
-- 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 exxpected 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 = { }
}
-- FiniteElementSpace1
FiniteElementSpace1 = {
-- Index of the god of dof into which the finite element space is defined.
......
......@@ -179,6 +179,39 @@ Domain3 = {
geometric_element_type_list = { }
}
Domain4 = {
-- 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 exxpected 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 = { }
}
-- FiniteElementSpace1
FiniteElementSpace1 = {
-- Index of the god of dof into which the finite element space is defined.
......
......@@ -46,7 +46,8 @@ namespace HappyHeart
{
highest_dimension = 1,
horizontal_source = 2,
spiral_source = 3
spiral_source = 3,
full_mesh = 4
};
......@@ -116,7 +117,8 @@ namespace HappyHeart
InputParameter::Domain<EnumUnderlyingType(DomainIndex::highest_dimension)>,
InputParameter::Domain<EnumUnderlyingType(DomainIndex::horizontal_source)>,
InputParameter::Domain<EnumUnderlyingType(DomainIndex::spiral_source)>,
InputParameter::Domain<EnumUnderlyingType(DomainIndex::full_mesh)>,
InputParameter::FEltSpace<EnumUnderlyingType(FEltSpaceIndex::highest_dimension)>,
InputParameter::FEltSpace<EnumUnderlyingType(FEltSpaceIndex::horizontal_source)>,
InputParameter::FEltSpace<EnumUnderlyingType(FEltSpaceIndex::spiral_source)>,
......
......@@ -18,6 +18,7 @@
# include "HappyHeart/Utilities/Filesystem/Folder.hpp"
# include "HappyHeart/Geometry/Domain/Domain.hpp"
# include "HappyHeart/Geometry/Domain/DomainManager.hpp"
# include "HappyHeart/FormulationSolver/VariationalFormulation.hpp"
......
......@@ -55,12 +55,14 @@ namespace HappyHeart
template <Advanced::ReactionLawNS::ReactionLawName ReactionLawNameT>
void ReactionDiffusionVariationalFormulation<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;
diffusion_tensor_ = InitScalarParameterFromInputData<Diffusion::Tensor<EnumUnderlyingType(TensorIndex::diffusion_tensor)>>("Diffusion tensor", mesh, input_parameter_data);
density_ = InitScalarParameterFromInputData<Diffusion::Density>("Density", mesh, input_parameter_data);
diffusion_tensor_ =
InitScalarParameterFromInputData<Diffusion::Tensor<EnumUnderlyingType(TensorIndex::diffusion_tensor)>>("Diffusion tensor",
domain,
input_parameter_data);
density_ = InitScalarParameterFromInputData<Diffusion::Density>("Density", domain, input_parameter_data);
if (!GetDiffusionTensor().IsConstant())
throw Exception("Current Reaction-Diffusion model is restricted to a constant diffusion tensor.", __FILE__, __LINE__);
......@@ -112,7 +114,9 @@ namespace HappyHeart
void ReactionDiffusionVariationalFormulation<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& felt_space_horizontal_source = god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::horizontal_source));
const auto& felt_space_spiral_source = god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::spiral_source));
......@@ -131,7 +135,7 @@ namespace HappyHeart
using parameter_type = InputParameter::TransientSource<EnumUnderlyingType(ForceIndexList::horizontal_source)>;
horizontal_source_parameter_ =
InitThreeDimensionalParameterFromInputData<parameter_type>("Horizontal source",
mesh,
domain,
input_parameter_data);
if (horizontal_source_parameter_ != nullptr)
......@@ -146,7 +150,7 @@ namespace HappyHeart
using parameter_type = InputParameter::TransientSource<EnumUnderlyingType(ForceIndexList::spiral_source)>;
spiral_source_parameter_ =
InitThreeDimensionalParameterFromInputData<parameter_type>("Spiral source",
mesh,
domain,
input_parameter_data);
if (spiral_source_parameter_ != nullptr)
......@@ -159,7 +163,7 @@ namespace HappyHeart
}
reaction_law_ = std::make_unique<reaction_law_type>(input_parameter_data,
mesh,
domain,
this->GetTimeManager(),
felt_space_highest_dimension.GetQuadratureRulePerTopology());
......
......@@ -82,7 +82,7 @@ namespace HappyHeart
*/
template <class InputParameterDataT>
explicit ReactionLaw(const InputParameterDataT& input_parameter_data,
const GeometricMeshRegion& mesh,
const Domain& domain,
const TimeManager& time_manager,
const QuadratureRulePerTopology& default_quadrature_rule_set);
......
......@@ -25,7 +25,7 @@ namespace HappyHeart
template <class InputParameterDataT>
ReactionLaw<ReactionLawName::CourtemancheRamirezNattel>::ReactionLaw(const InputParameterDataT& input_parameter_data,
const GeometricMeshRegion& mesh,
const Domain& domain,
const TimeManager& time_manager,
const QuadratureRulePerTopology& default_quadrature_rule_set)
: time_manager_(time_manager)
......@@ -64,213 +64,212 @@ namespace HappyHeart
constexpr const double value_at_quad_jsr = 1.49;
parameter_list_[EnumUnderlyingType(parameter_index::m)] =
std::make_unique<ScalarParameterAtQuadPt>("m",
mesh,
domain,
default_quadrature_rule_set,
value_at_quad_m,
this->GetTimeManager());
parameter_list_[EnumUnderlyingType(parameter_index::h)] =
std::make_unique<ScalarParameterAtQuadPt>("h",
mesh,
domain,
default_quadrature_rule_set,
value_at_quad_h,
this->GetTimeManager());
parameter_list_[EnumUnderlyingType(parameter_index::j)] =
std::make_unique<ScalarParameterAtQuadPt>("j",
mesh,
domain,
default_quadrature_rule_set,
value_at_quad_j,
this->GetTimeManager());
parameter_list_[EnumUnderlyingType(parameter_index::ao)] =
std::make_unique<ScalarParameterAtQuadPt>("ao",
mesh,
domain,
default_quadrature_rule_set,
value_at_quad_ao,
this->GetTimeManager());
parameter_list_[EnumUnderlyingType(parameter_index::io)] =
std::make_unique<ScalarParameterAtQuadPt>("io",
mesh,
domain,
default_quadrature_rule_set,
value_at_quad_io,
this->GetTimeManager());
parameter_list_[EnumUnderlyingType(parameter_index::ua)] =
std::make_unique<ScalarParameterAtQuadPt>("ua",
mesh,
domain,
default_quadrature_rule_set,
value_at_quad_ua,
this->GetTimeManager());
parameter_list_[EnumUnderlyingType(parameter_index::ui)] =
std::make_unique<ScalarParameterAtQuadPt>("ui",
mesh,
domain,
default_quadrature_rule_set,
value_at_quad_ui,
this->GetTimeManager());
parameter_list_[EnumUnderlyingType(parameter_index::xr)] =
std::make_unique<ScalarParameterAtQuadPt>("xr",
mesh,
domain,
default_quadrature_rule_set,
value_at_quad_xr,
this->GetTimeManager());
parameter_list_[EnumUnderlyingType(parameter_index::xs)] =
std::make_unique<ScalarParameterAtQuadPt>("xs",
mesh,
domain,
default_quadrature_rule_set,
value_at_quad_xs,
this->GetTimeManager());
parameter_list_[EnumUnderlyingType(parameter_index::d)] =
std::make_unique<ScalarParameterAtQuadPt>("d",
mesh,
domain,
default_quadrature_rule_set,
value_at_quad_d,
this->GetTimeManager());
parameter_list_[EnumUnderlyingType(parameter_index::f)] =
std::make_unique<ScalarParameterAtQuadPt>("f",
mesh,
domain,
default_quadrature_rule_set,
value_at_quad_f,
this->GetTimeManager());
parameter_list_[EnumUnderlyingType(parameter_index::fca)] =
std::make_unique<ScalarParameterAtQuadPt>("fca",
mesh,
domain,
default_quadrature_rule_set,
value_at_quad_fca,
this->GetTimeManager());
parameter_list_[EnumUnderlyingType(parameter_index::urel)] =
std::make_unique<ScalarParameterAtQuadPt>("urel",
mesh,
domain,
default_quadrature_rule_set,
value_at_quad_urel,
this->GetTimeManager());
parameter_list_[EnumUnderlyingType(parameter_index::vrel)] =
std::make_unique<ScalarParameterAtQuadPt>("vrel",
mesh,
domain,
default_quadrature_rule_set,
value_at_quad_vrel,
this->GetTimeManager());
parameter_list_[EnumUnderlyingType(parameter_index::wrel)] =
std::make_unique<ScalarParameterAtQuadPt>("wrel",
mesh,
domain,
default_quadrature_rule_set,
value_at_quad_wrel,
this->GetTimeManager());
parameter_list_[EnumUnderlyingType(parameter_index::nai)] =
std::make_unique<ScalarParameterAtQuadPt>("nai",
mesh,
domain,
default_quadrature_rule_set,
value_at_quad_nai,
this->GetTimeManager());
parameter_list_[EnumUnderlyingType(parameter_index::nao)] =
std::make_unique<ScalarParameterAtQuadPt>("nao",
mesh,
domain,
default_quadrature_rule_set,
value_at_quad_nao,
this->GetTimeManager());
parameter_list_[EnumUnderlyingType(parameter_index::cao)] =
std::make_unique<ScalarParameterAtQuadPt>("cao",
mesh,
domain,
default_quadrature_rule_set,
value_at_quad_cao,
this->GetTimeManager());
parameter_list_[EnumUnderlyingType(parameter_index::ki)] =
std::make_unique<ScalarParameterAtQuadPt>("ki",
mesh,
domain,
default_quadrature_rule_set,
value_at_quad_ki,
this->GetTimeManager());
parameter_list_[EnumUnderlyingType(parameter_index::ko)] =
std::make_unique<ScalarParameterAtQuadPt>("ko",
mesh,
domain,
default_quadrature_rule_set,
value_at_quad_ko,
this->GetTimeManager());
parameter_list_[EnumUnderlyingType(parameter_index::cai)] =
std::make_unique<ScalarParameterAtQuadPt>("cai",
mesh,
domain,
default_quadrature_rule_set,
value_at_quad_cai,
this->GetTimeManager());
parameter_list_[EnumUnderlyingType(parameter_index::naiont)] =
std::make_unique<ScalarParameterAtQuadPt>("naiont",
mesh,
domain,
default_quadrature_rule_set,
value_at_quad_naiont,
this->GetTimeManager());
parameter_list_[EnumUnderlyingType(parameter_index::kiont)] =
std::make_unique<ScalarParameterAtQuadPt>("kiont",
mesh,
domain,
default_quadrature_rule_set,
value_at_quad_kiont,
this->GetTimeManager());
parameter_list_[EnumUnderlyingType(parameter_index::caiont)] =
std::make_unique<ScalarParameterAtQuadPt>("caiont",
mesh,
domain,
default_quadrature_rule_set,
value_at_quad_caiont,
this->GetTimeManager());
parameter_list_[EnumUnderlyingType(parameter_index::ileak)] =
std::make_unique<ScalarParameterAtQuadPt>("ileak",
mesh,
domain,
default_quadrature_rule_set,
value_at_quad_ileak,
this->GetTimeManager());
parameter_list_[EnumUnderlyingType(parameter_index::iup)] =
std::make_unique<ScalarParameterAtQuadPt>("iup",
mesh,
domain,
default_quadrature_rule_set,
value_at_quad_iup,
this->GetTimeManager());
parameter_list_[EnumUnderlyingType(parameter_index::itr)] =
std::make_unique<ScalarParameterAtQuadPt>("itr",
mesh,
domain,
default_quadrature_rule_set,
value_at_quad_itr,
this->GetTimeManager());
parameter_list_[EnumUnderlyingType(parameter_index::irel)] =
std::make_unique<ScalarParameterAtQuadPt>("irel",
mesh,
domain,
default_quadrature_rule_set,
value_at_quad_irel,
this->GetTimeManager());
parameter_list_[EnumUnderlyingType(parameter_index::nsr)] =
std::make_unique<ScalarParameterAtQuadPt>("nsr",
mesh,
domain,
default_quadrature_rule_set,
value_at_quad_nsr,
this->GetTimeManager());
parameter_list_[EnumUnderlyingType(parameter_index::jsr)] =
std::make_unique<ScalarParameterAtQuadPt>("jsr",
mesh,
domain,
default_quadrature_rule_set,
value_at_quad_jsr,
this->GetTimeManager());
......
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