Commit 005cd69f authored by laura.pfeiffer's avatar laura.pfeiffer
Browse files

Merge remote-tracking branch 'origin/feature' into convergenceTestsMonodomain

parents 1621bfbe 3d4904fe
......@@ -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)
#---------------------------------------------------------------------------------------#
......
......@@ -9,10 +9,10 @@ h_min=1e-6
Model = ActiveStrainElasticity
MechModel = Static
MechDiscretization = Conforming
MechDiscretization = DG
CalculateReference = 1
ReferenceLevel = 4
CalculateReference = 0
ReferenceLevel = 2
InterpolateStartVector = false
MechProblem = LeftVentricle
......@@ -58,7 +58,7 @@ StretchFactorNormal = 4.0
######################
#ActiveMaterial = Linear
#ActiveMaterial = NeoHooke
#ActiveMaterial = Bonet
ActiveMaterial = Bonet
ActiveMaterial = Holzapfel
Density = 0.00
......@@ -77,8 +77,8 @@ PermeabilityPenalty = 20
#------------------------
# Linear Material
#------------------------
LinearMat_Lambda = 20
LinearMat_Mu = 10
LinearMat_Lambda = 0.6
LinearMat_Mu = 0.3
#------------------------
# NeoHooke Material
......@@ -89,10 +89,11 @@ NeoHooke_Mu = 10
#------------------------
# Guccione Material
#------------------------
GuccioneMat_C=10
GuccioneMat_bf=1
GuccioneMat_bs=1
GuccioneMat_bfs=1
GuccioneMat_C=0.313
GuccioneMat_bf=17.8
GuccioneMat_bs=7.1
GuccioneMat_bfs=12.4
#------------------------
......@@ -112,7 +113,7 @@ HolzapfelMat_bfs = 10
# Pressure Parameters#
######################
#ConstantPressure = -10
LVPressure = 10
LVPressure = 0.5 # 10
RVPressure = 0
RAPressure = 0
LAPressure = 0
......@@ -162,6 +163,9 @@ presmoothing = 3;
postsmoothing = 3;
DGSign=-1.0
DGPenalty= 25.0
Overlap = dG1
Overlap = NoOverlap
Overlap_Distribution = 0
......@@ -198,15 +202,15 @@ ElphyVerbose = -1
MechVerbose = -1
# = Plot Verbose =
PlotVTK = 0
PlotVTK = 1
PlottingSteps = 1
ElphySolverVTK = -1
PressureSolverVTK = -1
PressureSolverVTK = 1
DynamicSolverVTK = -1
CoupledSolverVTK = -1
PrestressPlotVerbose = 0
PressurePlotVerbose = -1
PressurePlotVerbose = 1
DynamicPlotVerbose = -1
......@@ -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);
......@@ -1246,7 +1246,51 @@ void DGElasticity::ProblemEvaluation(const Vector &u) const {
mout << "Apex Displacement = " << apexDisplacement << endl;
mout << "Apex Location = " << apexDisplacement[2] - 20.0 << endl;
mout << "Pint Displacement = " << pointDisplacemts << endl;
mout << "Point Displacement = " << pointDisplacemts << endl;
}
if(problemName =="LeftVentricle"){
VectorField apexDisplacement{};
VectorField pointDisplacemts{};
int cellCount{0};
int cellCount1{0};
for(cell c = u.cells(); c!= u.cells_end(); ++c){
for(int i = 0;i<c.Corners(); ++i){
if(c.Corner(i) == Point(84.434, -44.9479, -12.4509)){
DGVectorFieldElement E(u,c);
apexDisplacement += E.VectorValue(c.LocalCorner(i), u);
cellCount++;
}
}
}
cellCount = PPM->Sum(cellCount);
if(cellCount>0){
apexDisplacement = PPM->Sum(apexDisplacement) / cellCount;
}
const std::vector<Point> points = {Point(84.434, -44.9479, -12.4509),(-5.886750221,0,-16.16740036)};
for (row r = u.rows(); r != u.rows_end(); ++r) {
auto[i, isNear] = isIn(r(), points);
if (isNear) {
for (cell c = u.cells(); c != u.cells_end(); ++c) {
for (int i = 0; i < c.Corners(); ++i) {
if (c.Corner(i) == Point(84.434, -44.9479, -12.4509)) {
DGVectorFieldElement E(u, c);
pointDisplacemts += E.VectorValue(c.LocalCorner(i), u);
cellCount1++;
}
}
}
}
}
cellCount1 = PPM->Sum(cellCount1);
if(cellCount1>0){
pointDisplacemts = PPM->Sum(apexDisplacement) / cellCount;
}
mout << "Apex Displacement = " << apexDisplacement << endl;
mout << "Point Displacement = " << pointDisplacemts << endl;
}
}
......
......@@ -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;