Commit 2a78e345 authored by uheqb's avatar uheqb
Browse files

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

Merge branch '144-check-assemble-consistency' into 131-separate-finiteelasticity-from-lagrangeelasticity
parents 7c9303a0 44e08540
......@@ -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_TESTS ON)
set_option(BUILD_CARDIAC_TESTS ON)
set_option(BUILD_IBTLIBRARY OFF)
set_option(BUILD_TUTORIAL ON)
set_option(BUILD_TUTORIAL_TESTS ON)
#---------------------------------------------------------------------------------------#
......
......@@ -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,12 @@ void MainElasticity::Initialize() {
pressureProblem = getPressureProblem();
meshes = generateMesh();
initAssemble();
generateMatrixGraphs();
generateStartVectors();
initAssemble();
}
Vector &MainElasticity::Run(const Vector &start) {
......@@ -60,35 +61,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,22 +97,15 @@ 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);
}
int degree = 1;
bool isStatic = false;
config.get("MechPolynomialDegree", degree);
mechA = CreateFiniteElasticityAssemble(*elasticityProblem, *pressureProblem, *meshes, degree, isStatic);
auto prestress = false;
config.get("WithPrestress", prestress);
if (prestress) {
StaticPressureProblem prestressPressure(pressureProblem->InitialPressure());
LagrangeElasticity pEla(*elasticityProblem, prestressPressure);
LagrangeElasticity pEla(*elasticityProblem, prestressPressure, *meshes, degree,isStatic);
IterativePressureSolver prestressSolver(pEla, prestressPressure);
Vector p(0.0, *displacement);
prestressSolver.Solve(p, prestressSteps);
......
......@@ -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};
......
......@@ -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;
......
......@@ -3,6 +3,9 @@
#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) {
......@@ -121,3 +124,33 @@ 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);
if (elasticityModel == "DG") {
return std::make_unique<DGElasticity>(elasticityProblem, pressureProblem, mesh, degree, isStatic);
}
else if (elasticityModel == "EG") {
return std::make_unique<EGElasticity>(elasticityProblem, pressureProblem, mesh, degree, isStatic);
}
else {
return std::make_unique<LagrangeElasticity>(elasticityProblem, pressureProblem, 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);
}
}
......@@ -45,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;
......@@ -154,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;
}
......
......@@ -2,10 +2,17 @@
#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;
......
add_mpp_test(TestMaterial ELASTICITY)
add_mpp_test(TestCirculation ELASTICITY)
add_mpp_test(TestActiveStrain ELASTICITY)
add_mpp_test(TestAssembleConsistency ELASTICITY)
# These should not yet go in the pipeline
#add_mpp_test(TestLShape CardMech)
#add_mpp_test(TestDirichletBeam CardMech)
......
//
// Created by laura on 02.06.22.
//
#include "MainElasticity.hpp"
#include "TestEnvironment.hpp"
#include "TestConfigurations.hpp"
std::vector<std::string> elasticityProblems{"CardiacBeam"};
struct TestBeamParameter {
std::string discretization;
int problemLevel;
int degree;
std::string material;
std::string QuasiCompressiblePenalty;
double NewtonDamping;
std::string Overlap;
int OverlapDistribution;
int DGPenalty;
};
class ElasticityTest : public TestWithParam<TestBeamParameter> {
protected:
std::unique_ptr<ElasticityProblem> elastProblem;
std::unique_ptr<PressureProblem> pressProblem;
std::unique_ptr<FiniteElasticity> assemble;
std::unique_ptr<Meshes> meshes;
std::unique_ptr<Vector> u;
std::string model;
explicit ElasticityTest(const std::string &model) : model(model) {}
void SetUp() override {
std::map<string, string> testConfig = CARDIAC_TEST_CONFIG;
testConfig["Model"] = "ActiveStrainElasticity";
testConfig["MechModel"] = "Static";
testConfig["LAPressure"] = "0.004";
testConfig["Mesh"] = "BenchmarkBeam3DTet";
testConfig["MechProblem"] = "CardiacBeam"; //TODO noch anpassen auf mehrere Probleme
testConfig["VolumetricPenalty"] = "2";
testConfig["PermeabilityPenalty"] = "2";
testConfig["ConfigVerbose"] = "5";
testConfig["MechLevel"] = std::to_string(GetParam().problemLevel);
testConfig["MechDiscretization"] = GetParam().discretization;
testConfig["ActiveMaterial"] = GetParam().material;
testConfig["QuasiCompressiblePenalty"] = GetParam().QuasiCompressiblePenalty;
testConfig["NewtonDamping"] = std::to_string(GetParam().NewtonDamping);
testConfig["Overlap"] = GetParam().Overlap;
testConfig["Overlap_Distribution"] = std::to_string(GetParam().OverlapDistribution);
testConfig["DGPenalty"] = std::to_string(GetParam().DGPenalty);
/**************************************************
* Linear Material
**************************************************/
testConfig["LinearMat_Mu"] = "2";
testConfig["LinearMat_Lambda"] = "4";
/**************************************************
* Bonet Material
**************************************************/
testConfig["GuccioneMat_C"] = "2";
testConfig["GuccioneMat_bf"] = "8";
testConfig["GuccioneMat_bs"] = "2";
testConfig["GuccioneMat_bfs"] = "4";
Config::initialize(testConfig);
config.PrintInfo();
mout << "MechLevel Test " << GetParam().problemLevel << endl;
mout << "MechPolynomialDegree Test " << GetParam().degree << endl;
elastProblem = getElasticityProblem();//CardiacBeam
pressProblem = getPressureProblem();//"CardiacBeam"
meshes = MeshesCreator("BenchmarkBeam3DTet").CreateUnique();
assemble = CreateFiniteElasticityAssemble(*elastProblem, *pressProblem, model, *meshes,
GetParam().degree);
u = std::make_unique<Vector>(0.0, assemble->GetDisc());
assemble->Initialize(*u, 0.0, 1.0);
}
void TearDown() override {
Config::close();
}
};
/**************************************************
* Lagrange Consistency
**************************************************/
class LagrangeElasticityTest : public ElasticityTest {
protected:
LagrangeElasticityTest() : ElasticityTest(GetParam().discretization) {
}
};
TEST_P(LagrangeElasticityTest, AssembleConsistency) {
EXPECT_TRUE(isAssembleConsistent(*assemble, *u));
}
/*std::vector<TestBeamParameter> BeamParameters(int degree) {
return std::vector<TestBeamParameter>{
TestBeamParameter({"Conforming", degree}),
TestBeamParameter({"DG", degree}),
};
}*/
INSTANTIATE_TEST_SUITE_P(Degree_1_Consistency, LagrangeElasticityTest, Values(
TestBeamParameter({"Conforming", 0, 1, "Linear", "None", 1.0, "NoOverlap", 0})
/*TestBeamParameter({"Conforming", 1, 1, "Linear", "None", 1.0, "NoOverlap", 0 }),
TestBeamParameter({"Conforming", 2, 1, "Linear", "None", 1.0, "NoOverlap", 0}),
TestBeamParameter({"Conforming", 0, 1, "Bonet", "Ciarlet", 0.75, "NoOverlap", 0}),
TestBeamParameter({"Conforming", 1, 1, "Bonet", "Ciarlet", 0.75, "NoOverlap", 0}),
TestBeamParameter({"Conforming", 2, 1, "Bonet", "Ciarlet", 0.75, "NoOverlap", 0})*/
));
/*INSTANTIATE_TEST_SUITE_P(Degree_2_Consistency, LagrangeElasticityTest,Values(
TestBeamParameter({"Conforming", 0, 2, "Linear", "None", 1.0, "NoOverlap", 0}),
TestBeamParameter({"Conforming", 1, 2, "Linear", "None", 1.0, "NoOverlap", 0 }),
TestBeamParameter({"Conforming", 2, 2, "Linear", "None", 1.0, "NoOverlap", 0}),
TestBeamParameter({"Conforming", 0, 2, "Bonet", "Ciarlet", 0.75, "NoOverlap", 0}),
TestBeamParameter({"Conforming", 1, 2, "Bonet", "Ciarlet", 0.75, "NoOverlap", 0}),
TestBeamParameter({"Conforming", 2, 2, "Bonet", "Ciarlet", 0.75, "NoOverlap", 0})
));*/
/**************************************************
* DG Consistency
**************************************************/
class DGElasticityTest : public ElasticityTest {
protected:
DGElasticityTest() : ElasticityTest(GetParam().discretization) {
}
};
TEST_P(DGElasticityTest, AssembleConsistency) {
EXPECT_TRUE(isAssembleConsistent(*assemble, *u));
}
INSTANTIATE_TEST_SUITE_P(Degree_1_Consistency, DGElasticityTest,Values(
TestBeamParameter({"DG", 0, 1, "Linear", "None", 1.0, "dG1", 1, 9})
//TestBeamParameter({"DG", 1, 1, "Linear", "None", 1.0, "dG1", 1, 9 }),
//TestBeamParameter({"DG", 2, 1, "Linear", "None", 1.0, "dG1", 1, 9}),
//TestBeamParameter({"DG", 0, 1, "Bonet", "Ciarlet", 0.75, "dG1", 1, 9})
//TestBeamParameter({"DG", 1, 1, "Bonet", "Ciarlet", 0.75, "dG1", 1, 9})
//TestBeamParameter({"DG", 2, 1, "Bonet", "Ciarlet", 0.75, "dG1", 1, 9})
));
/*INSTANTIATE_TEST_SUITE_P(Degree_1_Consistency, DGElasticityTest,Values(
TestBeamParameter({"DG", 0, 2, "Linear", "None", 1.0, "dG1", 1, 9}),
//TestBeamParameter({"DG", 1, 2, "Linear", "None", 1.0, "dG1", 1, 9 }),
//TestBeamParameter({"DG", 2, 2, "Linear", "None", 1.0, "dG1", 1, 9}),
TestBeamParameter({"DG", 0, 2, "Bonet", "Ciarlet", 0.75, "dG1", 1, 9})
//TestBeamParameter({"DG", 1, 2, "Bonet", "Ciarlet", 0.75, "dG1", 1, 9})
//TestBeamParameter({"DG", 2, 2, "Bonet", "Ciarlet", 0.75, "dG1", 1, 9})
));*/
/**************************************************
* EG Consistency
**************************************************/
class EGElasticityTest : public ElasticityTest {
protected:
EGElasticityTest() : ElasticityTest(GetParam().discretization) {
}
};
TEST_P(EGElasticityTest, AssembleConsistency) {
EXPECT_TRUE(isAssembleConsistent(*assemble, *u));
}
INSTANTIATE_TEST_SUITE_P(Degree_1_Consistency, EGElasticityTest,Values(
TestBeamParameter({"EG", 0, 1, "Linear", "None", 1.0, "EG", 1, 2}),
//TestBeamParameter({"EG", 1, 1, "Linear", "None", 1.0, "dG1", 1, 2 }),
//TestBeamParameter({"EG", 2, 1, "Linear", "None", 1.0, "dG1", 1, 2}),
TestBeamParameter({"EG", 0, 1, "Bonet", "Ciarlet", 0.75, "EG", 1, 2})
//TestBeamParameter({"EG", 1, 1, "Bonet", "Ciarlet", 0.75, "dG1", 1, 2})
//TestBeamParameter({"EG", 2, 1, "Bonet", "Ciarlet", 0.75, "dG1", 1, 2})
));
/*INSTANTIATE_TEST_SUITE_P(Degree_1_Consistency, DGElasticityTest,Values(
TestBeamParameter({"EG", 0, 2, "Linear", "None", 1.0, "dG1", 1, 2}),
//TestBeamParameter({"EG", 1, 2, "Linear", "None", 1.0, "dG1", 1, 2 }),
//TestBeamParameter({"EG", 2, 2, "Linear", "None", 1.0, "dG1", 1, 2}),
TestBeamParameter({"EG", 0, 2, "Bonet", "Ciarlet", 0.75, "dG1", 1, 2})
//TestBeamParameter({"EG", 1, 2, "Bonet", "Ciarlet", 0.75, "dG1", 1, 2})
//TestBeamParameter({"EG", 2, 2, "Bonet", "Ciarlet", 0.75, "dG1", 1, 2})
));*/
int main(int argc, char **argv) {
return MppTest(
MppTestBuilder(argc, argv).
//WithSearchPath("/../../../..").
WithoutDefaultConfig().
WithScreenLogging().
WithPPM()
).RUN_ALL_MPP_TESTS();
}
\ No newline at end of file
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