Commit 7d1418b5 authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#1022 [STASH] Stuff before formatting turtle.

parent e828b8af
Pipeline #51476 failed with stage
in 0 seconds
......@@ -21,7 +21,7 @@ transient = {
-- Maximum time, if set to zero run a static case.
-- Expected format: VALUE
-- Constraint: v >= 0.
timeMax = .15,
timeMax = 0.,
} -- transient
......@@ -59,7 +59,7 @@ Mesh1 = {
-- Path of the mesh file to use.
-- Expected format: "VALUE"
mesh = "${HOME}/Codes/HappyHeart/Data/Mesh/elasticity_Nx50_Ny20_force_label.mesh",
mesh = "${HOME}/Codes/HappyHeart/Data/Mesh/elasticity_Nx3_Ny2_force_label.mesh",
-- Format of the input mesh.
-- Expected format: "VALUE"
......@@ -272,7 +272,7 @@ FiniteElementSpace1 = {
-- List of the shape function to use for each unknown;
-- Expected format: {"VALUE1", "VALUE2", ...}
shape_function_list = {"P2"},
shape_function_list = {"P1"},
-- List of the numbering subset to use for each unknown;
-- Expected format: {VALUE1, VALUE2, ...}
......@@ -297,7 +297,7 @@ FiniteElementSpace2 = {
-- List of the shape function to use for each unknown;
-- Expected format: {"VALUE1", "VALUE2", ...}
shape_function_list = {"P2"},
shape_function_list = {"P1"},
-- List of the numbering subset to use for each unknown;
-- Expected format: {VALUE1, VALUE2, ...}
......@@ -359,7 +359,7 @@ Solid = {
-- Value of the volumic mass in the case nature is 'constant'(and also initial value if nature is 'at_quadrature_point'; irrelevant otherwise).
-- Expected format: VALUE
-- Constraint: v > 0.
scalar_value = 1.,
scalar_value = 1.3,
-- Domain indices of the parameter in the case nature is 'piecewise_constant_by_domain'.
-- Expected format: {VALUE1, VALUE2, ...}
......@@ -387,12 +387,12 @@ Solid = {
-- How is given the Young modulus (as a constant, as a Lua function, per quadrature point, etc...)
-- Expected format: "VALUE"
-- Constraint: ops_in(v, {'constant', 'lua_function', 'piecewise_constant_by_domain'})
nature = "piecewise_constant_by_domain",
nature = "ignore",
-- Value of the Young modulus in the case nature is 'constant'(and also initial value if nature is 'at_quadrature_point'; irrelevant otherwise).
-- Expected format: VALUE
-- Constraint: v > 0.
scalar_value = 21e5,
scalar_value =8307692.023668617,
-- Domain indices of the parameter in the case nature is 'piecewise_constant_by_domain'.
-- Expected format: {VALUE1, VALUE2, ...}
......@@ -419,12 +419,12 @@ Solid = {
-- How is given the Poisson coefficient (as a constant, as a Lua function, per quadrature point, etc...)
-- Expected format: "VALUE"
-- Constraint: ops_in(v, {'constant', 'lua_function', 'piecewise_constant_by_domain'})
nature = "piecewise_constant_by_domain",
nature = "ignore",
-- Value of the Poisson coefficient in the case nature is 'constant'(and also initial value if nature is 'at_quadrature_point'; irrelevant otherwise).
-- Expected format: VALUE
-- Constraint: v > 0.
scalar_value = 0.3,
scalar_value = 0.03846150295857715,
-- Domain indices of the parameter in the case nature is 'piecewise_constant_by_domain'.
-- Expected format: {VALUE1, VALUE2, ...}
......@@ -453,12 +453,12 @@ Solid = {
-- How is given the kappa 1 (as a constant, as a Lua function, per quadrature point, etc...)
-- Expected format: "VALUE"
-- Constraint: ops_in(v, {'constant', 'lua_function', 'piecewise_constant_by_domain'})
nature = "piecewise_constant_by_domain",
nature = "constant",
-- Value of the kappa 1 in the case nature is 'constant'(and also initial value if nature is 'at_quadrature_point'; irrelevant otherwise).
-- Expected format: VALUE
-- Constraint: v > 0.
scalar_value = 500,
scalar_value = 0.5e3,
-- Domain indices of the parameter in the case nature is 'piecewise_constant_by_domain'.
-- Expected format: {VALUE1, VALUE2, ...}
......@@ -483,12 +483,12 @@ Solid = {
-- How is given the kappa 2 (as a constant, as a Lua function, per quadrature point, etc...)
-- Expected format: "VALUE"
-- Constraint: ops_in(v, {'constant', 'lua_function', 'piecewise_constant_by_domain'})
nature = "piecewise_constant_by_domain",
nature = "constant",
-- Value of the kappa 2 in the case nature is 'constant'(and also initial value if nature is 'at_quadrature_point'; irrelevant otherwise).
-- Expected format: VALUE
-- Constraint: v > 0.
scalar_value = 403346.1538461538,
scalar_value = 1.5e3,
-- Domain indices of the parameter in the case nature is 'piecewise_constant_by_domain'.
-- Expected format: {VALUE1, VALUE2, ...}
......@@ -513,12 +513,12 @@ Solid = {
-- How is given the parameter (as a constant, as a Lua function, per quadrature point, etc...)
-- Expected format: "VALUE"
-- Constraint: ops_in(v, {'ignore', 'constant', 'lua_function', 'at_quadrature_point', 'piecewise_constant_by_domain'})
nature = 'piecewise_constant_by_domain',
nature = 'constant',
-- Value of the volumic mass in the case nature is 'constant'(and also initial value if nature is
-- 'at_quadrature_point'; irrelevant otherwise).
-- Expected format: VALUE
scalar_value = 1750000.0,
scalar_value = 1750000,
-- Value of the volumic mass in the case nature is 'lua_function'(and also initial value if nature is
-- 'at_quadrature_point'; irrelevant otherwise).
......@@ -549,7 +549,7 @@ Solid = {
-- How is given the parameter (as a constant, as a Lua function, per quadrature point, etc...)
-- Expected format: "VALUE"
-- Constraint: ops_in(v, {'ignore', 'constant', 'lua_function', 'at_quadrature_point', 'piecewise_constant_by_domain'})
nature = 'constant',
nature = 'ignore',
-- Value of the volumic mass in the case nature is 'constant'(and also initial value if nature is
-- 'at_quadrature_point'; irrelevant otherwise).
......@@ -587,7 +587,7 @@ Solid = {
-- How is given the parameter (as a constant, as a Lua function, per quadrature point, etc...)
-- Expected format: "VALUE"
-- Constraint: ops_in(v, {'ignore', 'constant', 'lua_function', 'at_quadrature_point', 'piecewise_constant_by_domain'})
nature = 'constant',
nature = 'ignore',
-- Value of the volumic mass in the case nature is 'constant'(and also initial value if nature is
-- 'at_quadrature_point'; irrelevant otherwise).
......
......@@ -64,11 +64,11 @@
<CommandLineArguments>
<CommandLineArgument
argument = "--input_parameters ${HOME}/Codes/HappyHeart/Data/Lua/demo_input_hyperelasticity.lua"
isEnabled = "NO">
isEnabled = "YES">
</CommandLineArgument>
<CommandLineArgument
argument = "--input_parameters ${HOME}/Codes/HappyHeart/Data/Lua/demo_input_hyperelasticity_3d.lua"
isEnabled = "YES">
isEnabled = "NO">
</CommandLineArgument>
</CommandLineArguments>
<AdditionalOptions>
......
from ComparePoromechanics import *
if __name__ == "__main__":
freefem_output_dir = "/Users/sebastien/Codes/Freefem/Elasticity/Results"
hh_output_dir = "/Volumes/Data/sebastien/Sandbox/HYperelasticity"
for iteration in range(0, 3):
CompareMatrix("Tangent - {its}".format(its = iteration), \
"{dir}/tangent_{its}.m".format(dir = hh_output_dir, its = iteration + 1), \
"{dir}/tangent_{its}.txt".format(dir = freefem_output_dir, its = iteration), \
do_consider_matrices = True,
precision = 1.e-10,
do_padd_with_zero = True)
CompareVector("RHS - {its}".format(its = iteration), \
"{dir}/residual_{its}.m".format(dir = hh_output_dir, its = iteration), \
"{dir}/residual_{its}.txt".format(dir = freefem_output_dir, its = iteration), \
precision = 1.e-12,
do_padd_with_zero = False)
CompareVector("Solution - {its}".format(its = iteration), \
"{dir}/result_{its}.m".format(dir = hh_output_dir, its = iteration+1), \
"{dir}/result_{its}.txt".format(dir = freefem_output_dir, its = iteration), \
precision = 1.e-12,
do_padd_with_zero = False)
print("\n\n[WARNING] These tests are quick-and-dirty and might not account for discrepancies in very low values.")
from ComparePoromechanics import *
if __name__ == "__main__":
freefem_output_dir = "/Volumes/Data/sebastien/Sandbox/FreefemPoromechanics"
hh_output_dir = "/Volumes/Data/sebastien/HappyHeart/Results/Poromechanics"
CompareMatrix("T11", \
"{dir}/HH_T11.m".format(dir = hh_output_dir), \
"{dir}/NEW_T11.txt".format(dir = freefem_output_dir), \
do_consider_matrices = True,
precision = 1.e-12,
do_padd_with_zero = False)
CompareMatrix("T12", \
"{dir}/HH_T12.m".format(dir = hh_output_dir), \
"{dir}/NEW_T12.txt".format(dir = freefem_output_dir), \
do_consider_matrices = True,
precision = 1.e-12,
do_padd_with_zero = False)
CompareMatrix("T22", \
"{dir}/HH_MONO_T22.m".format(dir = hh_output_dir), \
"{dir}/MONO_T22.txt".format(dir = freefem_output_dir), \
do_consider_matrices = True,
precision = 1.e-10,
do_padd_with_zero = True)
CompareMatrix("T21", \
"{dir}/HH_T21.m".format(dir = hh_output_dir), \
"{dir}/NEW_T21.txt".format(dir = freefem_output_dir), \
do_consider_matrices = True,
precision = 1.e-10,
do_padd_with_zero = True)
CompareVector("Mixt variables", \
"{dir}/NEW_MIXT.m".format(dir = hh_output_dir), \
"{dir}/NEW_MIXT.txt".format(dir = freefem_output_dir), \
precision = 1.e-20,
do_padd_with_zero = True)
CompareMatrix("Monolithic matrix", \
"{dir}/NEW_MONOLITHIC.m".format(dir = hh_output_dir), \
"{dir}/NEW_MONOLITHIC.txt".format(dir = freefem_output_dir), \
do_consider_matrices = True,
precision = 1.e-10,
do_padd_with_zero = True)
CompareVector("RHS", \
"{dir}/NEW_RHS.m".format(dir = hh_output_dir), \
"{dir}/NEW_RHS.txt".format(dir = freefem_output_dir), \
precision = 1.e-12,
do_padd_with_zero = False)
CompareVector("FIRST IT SOL", \
"{dir}/NEW_SOL0.m".format(dir = hh_output_dir), \
"{dir}/NEW_SOL0.txt".format(dir = freefem_output_dir), \
precision = 1.e-12,
do_padd_with_zero = False)
print("\n\n[WARNING] These tests are quick-and-dirty and might not account for discrepancies in very low values.")
......@@ -378,7 +378,7 @@ namespace HappyHeart
const Unknown& GetSolidDisplacement() const noexcept;
private:
public:
//! Method that points to latest value of the sought quantity within the Newton loop.
const GlobalVector& GetSnesEvaluationState() const noexcept;
......
......@@ -230,25 +230,86 @@ namespace HappyHeart
}
template
<
class VariationalFormulationT
>
PetscErrorCode Viewer(SNES snes, PetscInt its, PetscReal norm, void* context_as_void)
{
static_cast<void>(snes);
assert(context_as_void);
const VariationalFormulationT* variational_formulation_ptr =
static_cast<const VariationalFormulationT*>(context_as_void);
assert(!(!variational_formulation_ptr));
const auto& variational_formulation = *variational_formulation_ptr;
const auto& vector = variational_formulation.GetSnesEvaluationState();
const PetscReal evaluation_state_min = vector.Min(__FILE__, __LINE__).second;
const PetscReal evaluation_state_max = vector.Max(__FILE__, __LINE__).second;
std::ostringstream oconv;
oconv << VariationalFormulationT::ClassName() << " --- Norm is %14.12e and extrema are %8.6e and %8.6e\n";
Wrappers::Petsc::PrintMessageOnFirstProcessor(oconv.str().c_str(),
variational_formulation.MpiHappyHeart(),
__FILE__, __LINE__,
its,
norm,
evaluation_state_min,
evaluation_state_max
);
decltype(auto) numbering_subset = variational_formulation.GetNumberingSubset();
decltype(auto) tangent = variational_formulation.GetSystemMatrix(numbering_subset, numbering_subset);
decltype(auto) residual = variational_formulation.GetSystemRhs(numbering_subset);
decltype(auto) sol = variational_formulation.GetSystemSolution(numbering_subset);
tangent.View(variational_formulation.MpiHappyHeart(),
"/Volumes/Data/sebastien/Sandbox/HYperelasticity/tangent_" + std::to_string(its) + ".m",
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
residual.View(variational_formulation.MpiHappyHeart(),
"/Volumes/Data/sebastien/Sandbox/HYperelasticity/residual_" + std::to_string(its) + ".m",
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
sol.View(variational_formulation.MpiHappyHeart(),
"/Volumes/Data/sebastien/Sandbox/HYperelasticity/result_" + std::to_string(its) + ".m",
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
return 0;
}
template
<
class HyperelasticLawT,
HyperelasticityNS::TimeScheme TimeSchemeT
>
inline Wrappers::Petsc::Snes::SNESViewer VariationalFormulationHyperElasticity<HyperelasticLawT, TimeSchemeT>
::ImplementSnesViewer() const
inline Wrappers::Petsc::Snes::SNESViewer
VariationalFormulationHyperElasticity<HyperelasticLawT, TimeSchemeT>::ImplementSnesViewer() const
{
return SolverNS::SnesInterface<self>::Viewer;
return Viewer<VariationalFormulationHyperElasticity<HyperelasticLawT, TimeSchemeT>>;
}
template
<
class HyperelasticLawT,
HyperelasticityNS::TimeScheme TimeSchemeT
>
inline Wrappers::Petsc::Snes::SNESConvergenceTestFunction
VariationalFormulationHyperElasticity<HyperelasticLawT, TimeSchemeT>::ImplementSnesConvergenceTestFunction() const
VariationalFormulationHyperElasticity<HyperelasticLawT, TimeSchemeT>
::ImplementSnesConvergenceTestFunction() const
{
return nullptr;
}
......@@ -350,8 +411,9 @@ namespace HappyHeart
break;
}
this->template ApplyEssentialBoundaryCondition<VariationalFormulationNS::On::system_rhs>(numbering_subset,
numbering_subset);
this->template ApplyEssentialBoundaryCondition<VariationalFormulationNS::On::system_rhs,
BoundaryConditionMethod::penalization>(numbering_subset,
numbering_subset);
}
......@@ -376,8 +438,9 @@ namespace HappyHeart
break;
}
this->template ApplyEssentialBoundaryCondition<VariationalFormulationNS::On::system_matrix>(numbering_subset,
numbering_subset);
this->template ApplyEssentialBoundaryCondition<VariationalFormulationNS::On::system_matrix,
BoundaryConditionMethod::penalization>(numbering_subset,
numbering_subset);
}
......
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