Commit c0f1b7d4 authored by laura.stengel's avatar laura.stengel
Browse files

Merge branch '144-check-assemble-consistency' into 'feature'

Resolve "Check Assemble Consistency"

Closes #144

See merge request !160
parents 81e45e4f 21613a2b
Pipeline #219308 passed with stages
in 32 minutes and 10 seconds
......@@ -5,16 +5,20 @@ include(mpp/CMakeOptionMacros.txt)
add_compile_options(-Werror=return-type)
set_option(COMPILER_OPTIMIZE -O3)
#set_option(COMPILER_OPTIMIZE -O3)
set_option(SPACE_DIM 3)
set_option(USE_DATAMESH ON)
set_option(MPP_BUILD_TYPE MppRelease)
set_option(MPP_BUILD_TYPE MppDebug)
#set_option(MPP_BUILD_TYPE MppRelease)
set_option(SUPPRESS_WARNINGS OFF)
set_option(AFFINE_LINEAR_TRAFO ON)
set_option(BUILD_TESTS OFF)
set_option(BUILD_CARDIAC_TESTS ON)
set_option(BUILD_IBTLIBRARY OFF)
set_option(BUILD_TUTORIAL OFF)
set_option(BUILD_TUTORIAL_TESTS OFF)
#---------------------------------------------------------------------------------------#
......
......@@ -189,16 +189,16 @@ void MainCoupled::initAssemble() {
auto elasticityModel = elasticityProblem->Model();
if (elasticityModel == "DG") {
mechA = std::make_unique<DGElasticity>(*elasticityProblem, *pressureProblem);
mechA = CreateFiniteElasticityAssemble(*elasticityProblem, *pressureProblem, *mechMeshes, mechDegree);
} else {
mechA = std::make_unique<LagrangeElasticity>(*elasticityProblem, *pressureProblem);
mechA = CreateFiniteElasticityAssemble(*elasticityProblem, *pressureProblem, *mechMeshes, mechDegree);
}
auto prestress = false;
config.get("WithPrestress", prestress);
if (prestress) {
StaticPressureProblem prestressPressure(pressureProblem->InitialPressure());
LagrangeElasticity pEla(*elasticityProblem, prestressPressure, true);
LagrangeElasticity pEla(*elasticityProblem, prestressPressure, *mechMeshes, mechDegree, true);
IterativePressureSolver prestressSolver(pEla, prestressPressure);
Vector p(0.0, *displacement);
prestressSolver.Solve(p, prestressSteps);
......
......@@ -19,11 +19,29 @@ void MainElasticity::Initialize() {
pressureProblem = getPressureProblem();
meshes = generateMesh();
initAssemble();
generateMatrixGraphs();
generateStartVectors();
initAssemble();
auto prestress = false;
config.get("WithPrestress", prestress);
if (prestress) {
int degree = 1;
bool isStatic = false;
config.get("MechPolynomialDegree", degree);
auto elasticityModel = elasticityProblem->Model();
StaticPressureProblem prestressPressure(pressureProblem->InitialPressure());
prestressedEla = CreateFiniteElasticityAssemble(*elasticityProblem, prestressPressure,*meshes, degree, isStatic);
IterativePressureSolver prestressSolver(*prestressedEla, prestressPressure);
Vector p(0.0, *displacement);
prestressSolver.Solve(p, prestressSteps);
mechA->InitializePrestress(p);
}
}
Vector &MainElasticity::Run(const Vector &start) {
......@@ -60,35 +78,12 @@ std::unique_ptr<Meshes> MainElasticity::generateMesh() const {
CreateUnique();
}
std::unique_ptr<IDiscretization> MainElasticity::generateDisc(const Meshes &mesh) const {
auto elasticityModel = elasticityProblem->Model();
if (elasticityModel == "DG") {
return std::make_unique<DGDiscretization>(*meshes, discDegree, 3);
}
if (elasticityModel == "dg") {
return std::make_unique<DGDiscretization>(*meshes, discDegree, 3);
}
if (elasticityModel == "EG") {
return std::make_unique<EGDiscretization>(*meshes, discDegree, 3);
}
return std::make_unique<LagrangeDiscretization>(mesh, discDegree, 3);
}
void MainElasticity::generateMatrixGraphs() {
auto elasticityModel = elasticityProblem->Model();
vectorDisc = generateDisc(*meshes);
if (elasticityModel == "DG") {
scalarDisc = std::make_unique<DGDiscretization>(*meshes, discDegree, 1);
dataDiscretization = std::make_unique<DGDiscretization>(*meshes, discDegree, 3);
}
else if (elasticityModel == "EG") {
scalarDisc = std::make_unique<EGDiscretization>(*meshes, discDegree, 1);
dataDiscretization = std::make_unique<EGDiscretization>(*meshes, discDegree, 3);
}
else{
scalarDisc = std::make_unique<LagrangeDiscretization>(*meshes, discDegree, 1);
dataDiscretization = std::make_unique<LagrangeDiscretization>(*meshes, discDegree, 3);
}
vectorDisc = mechA->CreateDisc(3);
scalarDisc = mechA->CreateDisc(1);
dataDiscretization = mechA->CreateDisc(3);
}
void MainElasticity::generateStartVectors() {
......@@ -119,27 +114,10 @@ Vector &MainElasticity::GetDisplacement() const {
}
void MainElasticity::initAssemble() {
auto elasticityModel = elasticityProblem->Model();
if (elasticityModel == "DG") {
mechA = std::make_unique<DGElasticity>(*elasticityProblem, *pressureProblem);
}
else if (elasticityModel == "EG") {
mechA = std::make_unique<EGElasticity>(*elasticityProblem, *pressureProblem);
}
else {
mechA = std::make_unique<LagrangeElasticity>(*elasticityProblem, *pressureProblem);
}
auto prestress = false;
config.get("WithPrestress", prestress);
if (prestress) {
StaticPressureProblem prestressPressure(pressureProblem->InitialPressure());
LagrangeElasticity pEla(*elasticityProblem, prestressPressure);
IterativePressureSolver prestressSolver(pEla, prestressPressure);
Vector p(0.0, *displacement);
prestressSolver.Solve(p, prestressSteps);
mechA->InitializePrestress(p);
}
int degree = 1;
bool isStatic = false;
config.get("MechPolynomialDegree", degree);
mechA = std::make_unique<LagrangeElasticity>(*elasticityProblem, *pressureProblem, *meshes, degree, isStatic);
}
double MainElasticity::EvaluateQuantity(const Vector &solution, const std::string &quantity) const {
......
......@@ -12,7 +12,6 @@ class MainElasticity : public MainCardMech {
std::unique_ptr<Meshes> generateMesh() const;
std::unique_ptr<IDiscretization> generateDisc(const Meshes &mesh) const;
int discDegree{1};
......@@ -24,6 +23,7 @@ class MainElasticity : public MainCardMech {
std::unique_ptr<IDiscretization> scalarDisc = nullptr;
std::unique_ptr<FiniteElasticity> mechA = nullptr;
std::unique_ptr<FiniteElasticity> prestressedEla = nullptr;
std::unique_ptr<TimeSeries> tSeries = nullptr;
......
......@@ -5,8 +5,8 @@
#include "DGVectorFieldFaceElement.hpp"
#include "VectorAccess.hpp"
DGElasticity::DGElasticity(const ElasticityProblem &eP, const PressureProblem &pP)
: FiniteElasticity(eP, pP) {
DGElasticity::DGElasticity(const ElasticityProblem &eP, const PressureProblem &pP, const Meshes &mesh, int degree, bool isStatic)
: FiniteElasticity(eP, pP), disc(mesh, degree, size) {
config.get("DGSign", sign);
config.get("DGPenalty", penalty);
config.get("MechPolynomialDegree", degree);
......
......@@ -3,18 +3,26 @@
#include "FiniteElasticity.hpp"
#include "DGDiscretization.hpp"
class DGElasticity : public FiniteElasticity {
double sign{1};
double sign{-1};
double penalty{1};
double degree{1};
int size{3};
const DGDiscretization disc;
void dirichletBnd(const cell &c, Vector &U, int face) const;
void fixationBnd(const cell &c, Vector &U, int face) const;
public:
DGElasticity(const ElasticityProblem &eP, const PressureProblem &pP);
DGElasticity(const ElasticityProblem &eP, const PressureProblem &pP, const Meshes &mesh, int degree, bool isStatic = false);
std::unique_ptr<IDiscretization> CreateDisc(int sizes) const override {return std::make_unique<DGDiscretization>(disc.GetMeshes(), degree,sizes);}
const IDiscretization &GetDisc() const override { return disc; }
void Initialize(Vector &u) const override;
......
......@@ -13,8 +13,8 @@
#include "EGVectorFieldElement.hpp"
#include "EGVectorFieldFaceElement.hpp"
EGElasticity::EGElasticity(const ElasticityProblem &eP, const PressureProblem &pP)
: FiniteElasticity(eP, pP) {
EGElasticity::EGElasticity(const ElasticityProblem &eP, const PressureProblem &pP, const Meshes &mesh, int degree, bool isStatic)
: FiniteElasticity(eP, pP) , disc(mesh, degree, size){
config.get("DGSign", sign);
config.get("DGPenalty", penalty);
}
......
......@@ -2,17 +2,25 @@
#define EGELASTICITY_H
#include "FiniteElasticity.hpp"
#include "EGDiscretization.hpp"
class EGElasticity : public FiniteElasticity {
double sign{1};
double sign{-1};
double penalty{1};
int size{3};
const EGDiscretization disc;
void dirichletBnd(const cell &c, Vector &U, int face) const;
void fixationBnd(const cell &c, Vector &U, int face) const;
public:
EGElasticity(const ElasticityProblem &eP, const PressureProblem &pP);
EGElasticity(const ElasticityProblem &eP, const PressureProblem &pP , const Meshes &mesh, int degree, bool isStatic = false);
std::unique_ptr<IDiscretization> CreateDisc(int sizes) const override {return std::make_unique<EGDiscretization>(disc, degree,sizes);}
const IDiscretization &GetDisc() const override { return disc; }
void Initialize(Vector &u) const override;
......
#include <VectorFieldFaceElement.hpp>
#include <LagrangeDiscretization.hpp>
//#include <LagrangeDiscretization.hpp>
#include "Plotting.hpp"
#include "FiniteElasticity.hpp"
#include "VectorFieldElement.hpp"
#include "LagrangeElasticity.hpp"
#include "DGElasticity.hpp"
#include "EGElasticity.hpp"
FiniteElasticity::FiniteElasticity(const ElasticityProblem &eP, const PressureProblem &pP) :
CardiacAssemble(), eProblem(eP), pProblem(pP), check_p(false) {
......@@ -32,41 +35,6 @@ FiniteElasticity::FiniteElasticity(const ElasticityProblem &eP, const PressurePr
}
void
FiniteElasticity::dirichletBnd(const cell &c, const Vector &U, RowBndValues &u_c, int face) const {
ScalarElement elem(U, c);
for (int j = 0; j < U.NodalPointsOnFace(c, face); ++j) {
int k = U.NodalPointOnFace(c, face, j);
VectorField D(eProblem.Deformation(timeStep, elem.NodalPoint(k)));
for (int d = 0; d < c.dim(); ++d) {
u_c(k, d) = D[d];
u_c.D(k, d) = true;
}
}
}
void FiniteElasticity::planeFixationBnd(const cell &c, const Vector &U, RowBndValues &u_c,
int face, int d) const {
for (int j = 0; j < U.NodalPointsOnFace(c, face); ++j) {
int k = U.NodalPointOnFace(c, face, j);
u_c(k, d) = 0.0;
u_c.D(k, d) = true;
}
}
void
FiniteElasticity::fixationBnd(const cell &c, const Vector &U, RowBndValues &u_c,
int face) const {
for (int j = 0; j < U.NodalPointsOnFace(c, face); ++j) {
int k = U.NodalPointOnFace(c, face, j);
for (int d = 0; d < c.dim(); ++d) {
u_c(k, d) = 0.0;
u_c.D(k, d) = true;
}
}
}
void FiniteElasticity::Initialize(Vector &u, double t, double dt) {
currentTime = t;
timeStep = dt;
......@@ -123,8 +91,6 @@ void FiniteElasticity::Initialize(const cell &c, Vector &U) const {
}
}
void FiniteElasticity::Finalize(Vector &u, double t, double dt) {
UpdateAcceleration(u, dt);
}
......@@ -151,15 +117,6 @@ void FiniteElasticity::Initialize(Vectors &vecs, double t, double dt) {
UpdateStretch(vecs[1]);
}
void FiniteElasticity::UpdateStretch(const Vector &gamma_f) {
for (row r = gamma->rows(); r != gamma->rows_end(); ++r) {
double g = gamma_f(r, 0);
(*gamma)(r, 0) = g;
(*gamma)(r, 2) = factor_gn * g;
// To ensure det(FA)=1 [See Bonet et al.(2019)]
(*gamma)(r, 1) = (1.0 / (1.0 + g)) * (1.0 / (1.0 + factor_gn * g)) - 1.0;
}
}
void FiniteElasticity::Finalize(Vectors &vecs, double t, double dt) {
Finalize(vecs[0], t, dt);
......@@ -167,3 +124,24 @@ void FiniteElasticity::Finalize(Vectors &vecs, double t, double dt) {
GetInvariant(vecs[0], vecs[2]);
}
std::unique_ptr<FiniteElasticity> CreateFiniteElasticityAssemble(
ElasticityProblem &elasticityProblem, PressureProblem &pressureProblem, const Meshes &mesh, int degree, bool isStatic){
auto elasticityModel = elasticityProblem.Model();
if (elasticityModel.empty()) config.get("MechDiscretization", elasticityModel);
return CreateFiniteElasticityAssemble(elasticityProblem, pressureProblem, elasticityModel, mesh, degree, isStatic);
}
std::unique_ptr<FiniteElasticity> CreateFiniteElasticityAssemble(
ElasticityProblem &elasticityProblem, PressureProblem &pressureProblem, std::string modelName, const Meshes &mesh, int degree, bool isStatic){
if (modelName.empty()) config.get("MechDiscretization", modelName);
if (modelName == "DG") {
return std::make_unique<DGElasticity>(elasticityProblem, pressureProblem, mesh, degree, isStatic);
}
else if (modelName == "EG") {
return std::make_unique<EGElasticity>(elasticityProblem, pressureProblem, mesh, degree, isStatic);
}
else {
return std::make_unique<LagrangeElasticity>(elasticityProblem, pressureProblem, mesh, degree, isStatic);
}
}
......@@ -8,11 +8,17 @@
class FiniteElasticity : public CardiacAssemble {
void fixationBnd(const cell &c, const Vector &U, RowBndValues &u_c, int face) const;
virtual void fixationBnd(const cell &c, const Vector &U, RowBndValues &u_c, int face) const{
Warning("dirichletBnd is not defined")
};
void planeFixationBnd(const cell &c, const Vector &U, RowBndValues &u_c, int face, int d) const;
virtual void planeFixationBnd(const cell &c, const Vector &U, RowBndValues &u_c, int face, int d) const{
Warning("dirichletBnd is not defined")
};
void dirichletBnd(const cell &c, const Vector &U, RowBndValues &u_c, int face) const;
virtual void dirichletBnd(const cell &c, const Vector &U, RowBndValues &u_c, int face) const{
Warning("dirichletBnd is not defined")
};
protected:
const ElasticityProblem &eProblem;
......@@ -39,9 +45,18 @@ protected:
double factor_gs{1.0};
double factor_gn{1.0};
//Generate Discretization
int degree{1};
std::unique_ptr<Meshes> meshes = nullptr;
std::unique_ptr<ElasticityProblem> elasticityProblem = nullptr;
public:
FiniteElasticity(const ElasticityProblem &eP, const PressureProblem &pP);
virtual const IDiscretization &GetDisc() const = 0;
virtual std::unique_ptr<IDiscretization> CreateDisc(int size) const = 0;
const char *Name() const override { return "Cardiac Elasticity"; }
void Initialize(Vector &u, double t, double dt) override;
......@@ -124,7 +139,9 @@ public:
return 0.0;
}
void UpdateStretch(const Vector &gamma_f);
virtual void UpdateStretch(const Vector &gamma_f) const {
Warning("UpdateStretch is not defined")
};
virtual void SetDisplacement(Vector &u) const {
Warning("SetDisplacement is not defined")
......@@ -146,4 +163,12 @@ public:
}
};
std::unique_ptr<FiniteElasticity> CreateFiniteElasticityAssemble(
ElasticityProblem &elasticityProblem, PressureProblem &pressureProblem, const Meshes &mesh, int degree, bool isStatic = false);
std::unique_ptr<FiniteElasticity> CreateFiniteElasticityAssemble(
ElasticityProblem &elasticityProblem, PressureProblem &pressureProblem, std::string modelName, const Meshes &mesh, int degree, bool isStatic = false);
#endif //FINITEELASTICITY_H
......@@ -5,9 +5,9 @@
#include "LagrangeElasticity.hpp"
#include "LagrangeTransfer.hpp"
LagrangeElasticity::LagrangeElasticity(const ElasticityProblem &eP, const PressureProblem &pP,
LagrangeElasticity::LagrangeElasticity(const ElasticityProblem &eP, const PressureProblem &pP, const Meshes &mesh, int degree,
bool isStatic)
: FiniteElasticity(eP, pP) {
: FiniteElasticity(eP, pP) , disc(mesh, degree, size){
rho = isStatic ? 0.0 : rho;
}
......@@ -465,6 +465,39 @@ double LagrangeElasticity::EnergyError(const Vector &u, const Vector &reference)
return EnergyNorm(u - projection);
}
void LagrangeElasticity::planeFixationBnd(const cell &c, const Vector &U, RowBndValues &u_c,
int face, int d) const {
for (int j = 0; j < U.NodalPointsOnFace(c, face); ++j) {
int k = U.NodalPointOnFace(c, face, j);
u_c(k, d) = 0.0;
u_c.D(k, d) = true;
}
}
void
LagrangeElasticity::fixationBnd(const cell &c, const Vector &U, RowBndValues &u_c,
int face) const {
for (int j = 0; j < U.NodalPointsOnFace(c, face); ++j) {
int k = U.NodalPointOnFace(c, face, j);
for (int d = 0; d < c.dim(); ++d) {
u_c(k, d) = 0.0;
u_c.D(k, d) = true;
}
}
}
void LagrangeElasticity::dirichletBnd(const cell &c, const Vector &U, RowBndValues &u_c, int face) const {
ScalarElement elem(U, c);
for (int j = 0; j < U.NodalPointsOnFace(c, face); ++j) {
int k = U.NodalPointOnFace(c, face, j);
VectorField D(eProblem.Deformation(timeStep, elem.NodalPoint(k)));
for (int d = 0; d < c.dim(); ++d) {
u_c(k, d) = D[d];
u_c.D(k, d) = true;
}
}
}
void LagrangeElasticity::InitializePrestress(const Vector &U) {
Prestress = std::make_unique<Vector>(0.0, U);
Prestress->Clear();
......@@ -824,6 +857,16 @@ void LagrangeElasticity::SetDisplacement(Vector &u) const {
}
}
void LagrangeElasticity::UpdateStretch(const Vector &gamma_f) const {
for (row r = gamma->rows(); r != gamma->rows_end(); ++r) {
double g = gamma_f(r, 0);
(*gamma)(r, 0) = g;
(*gamma)(r, 2) = factor_gn * g;
// To ensure det(FA)=1 [See Bonet et al.(2019)]
(*gamma)(r, 1) = (1.0 / (1.0 + g)) * (1.0 / (1.0 + factor_gn * g)) - 1.0;
}
}
void LagrangeElasticity::InitializeTimeStepping(double t, double dt, const Vector &U) {
auto u = new Vector(U.GetDisc(), U.SpaceLevel(), U.TimeLevel());
u->Clear();
......
......@@ -2,10 +2,23 @@
#define CARDMECH_LAGRANGEELASTICITY_HPP
#include "FiniteElasticity.hpp"
#include "LagrangeDiscretization.hpp"
class LagrangeElasticity : public FiniteElasticity {
int size{3};
const LagrangeDiscretization disc;
public:
LagrangeElasticity(const ElasticityProblem &eP, const PressureProblem &pP, bool isStatic = false);
LagrangeElasticity(const ElasticityProblem &eP, const PressureProblem &pP, const Meshes &mesh, int degree, bool isStatic = false);
const IDiscretization &GetDisc() const override { return disc; }
std::unique_ptr<IDiscretization> CreateDisc(int sizes) const override {return std::make_unique<LagrangeDiscretization>(disc, degree,sizes);}
void fixationBnd(const cell &c, const Vector &U, RowBndValues &u_c, int face) const override;
void planeFixationBnd(const cell &c, const Vector &U, RowBndValues &u_c, int face, int d) const override;
void dirichletBnd(const cell &c, const Vector &U, RowBndValues &u_c, int face) const override;
void InitializePrestress(const Vector &u) override;
......@@ -96,6 +109,8 @@ public:
void SetDisplacement(Vector &u) const override;
void UpdateStretch(const Vector &gamma_f) const override;
VectorField CalculateStretch(double t, const Point &xi) const;
};
......
//
// Created by laura on 13.06.22.
//
#ifndef CARDMECH_TESTCARDIACBEAMCONFIGURATION_H
#define CARDMECH_TESTCARDIACBEAMCONFIGURATION_H
#include <map>
struct ConfigMaps {
std::string confPath;
std::map<std::string, std::string> defaultMap;
std::map<std::string, std::string> additionalMap;
};
static const std::map<std::string, std::string> CARDIACBEAM_TEST_CONFIG{
// ----- Problem Settings -----
{"GeoPath", "../../../../geo/"},
{"EulerReAssemble", "0"},
{"PressureIterations", "10"},
{"PressureDepth", "0"},
// === Linear Solver === //
{"LinearSolver", "gmres"},
{"LinearEpsilon", "1e-8"},
{"LinearReduction", "1e-8"},
{"LinearReduction", "1e-11"},
{"LinearSteps", "10000"},
{"BasePreconditioner", "GaussSeidel"},
{"Preconditioner", "GaussSeidel"},
// =================================================================
// === Gating Solver ===============================================
{"GatingSolver", "gmres"},
{"GatingPreconditioner", "GaussSeidel"},
{"GatingEpsilon", "1e-8"},
{"GatingReduction", "1e-12"},
// === Elphy Solver ================================================
{"ElphySolver", "gmres"},
{"ElphyPreconditioner", "GaussSeidel"},
{"ElphyEpsilon", "1e-8"},
{"ElphyReduction", "1e-12"},
// === Stretch Solver ==============================================
{"StretchSolver", "gmres"},
{"StretchPreconditioner", "GaussSeidel"},
{"StretchEpsilon", "1e-8"},
{"StretchReduction", "1e-12"},
// === Mechanics Solver ============================================
{"MechSolver", "gmres"},
{"MechPreconditioner", "GaussSeidel"},