Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • kit/mpp/cardmech
1 result
Show changes
Commits on Source (10)
Hide whitespace changes
Inline Side-by-side
Showing
with 415 additions and 36 deletions
mpp @ e448af94
Compare cfd2f66c...e448af94
Subproject commit cfd2f66c53e10c51415ddc429e56a002ad7b97b0 Subproject commit e448af94f8abb4d860b9ce064dd363c663353690
src/elasticity/MainElasticity.cpp
View file @ 8f9a221d
...@@ -96,9 +96,11 @@ Vector &MainElasticity::Run() { ...@@ -96,9 +96,11 @@ Vector &MainElasticity::Run() {
CalculatePrestress(*mechA, *displacement); CalculatePrestress(*mechA, *displacement);
} }
if (dyn == "Static") { if (dyn == "Static") {
mout << "=== Running MainElasticity after Prestress: ===============" << endl;
elasticityProblem->EvaluationResults(*displacement); elasticityProblem->EvaluationResults(*displacement);
// Generate IterativePressureSolver in src/elasticity/solvers/IterativePressureSolver.hpp // Generate IterativePressureSolver in src/elasticity/solvers/IterativePressureSolver.hpp
IterativePressureSolver mSolver(*elasticityProblem, pressureSteps); IterativePressureSolver mSolver(*elasticityProblem, pressureSteps);
*displacement = 0;
mSolver.Method(*mechA, *displacement); mSolver.Method(*mechA, *displacement);
} else { } else {
// Generate ElastodynamicTimeIntegrator in src/elasticity/solvers/DynamicMechanicsSolver.hpp // Generate ElastodynamicTimeIntegrator in src/elasticity/solvers/DynamicMechanicsSolver.hpp
......
src/elasticity/assemble/DGElasticity.cpp
View file @ 8f9a221d
...@@ -568,6 +568,44 @@ void DGElasticity::Jacobi(const cell &c, const Vector &U, Matrix &A) const { ...@@ -568,6 +568,44 @@ void DGElasticity::Jacobi(const cell &c, const Vector &U, Matrix &A) const {
} }
} }
//TODO stimmt das so?
void DGElasticity::TractionStiffness(Matrix &matrix) const {
matrix = 0;
for (cell c = matrix.cells(); c != matrix.cells_end(); ++c) {
DGVectorFieldElement E(matrix, *c);
DGSizedRowEntries A_c(matrix, *c, E.shape_size(), *c, E.shape_size());
for (int f = 0; f < c.Faces(); ++f) {
DGVectorFieldFaceElement faceElem(matrix, *c, f); //passt das hier zb?
if (matrix.OnBoundary(*c, f)) {
int bc = matrix.GetMesh().BoundaryFacePart(c.Face(f));
if ((bc >= ROBIN_BC_EPI) && (bc <= ROBIN_BC_BASE)) {
for (int q = 0; q < faceElem.nQ(); ++q) {
VectorField n = faceElem.QNormal(q);
auto N = Product(n, n);
double w = faceElem.QWeight(q);
for (int i = 0; i < faceElem.shape_size(); ++i) {
for (int k = 0; k < E.Dim(); ++k) {
auto phi_ik = faceElem.VectorComponentValue(q, i, k);
auto Nphi_ik = N * faceElem.VectorComponentValue(q, i, k);
for (int j = 0; j < faceElem.shape_size(); ++j) {
for (int l = 0; l < E.Dim(); ++l) {
auto phi_jl = faceElem.VectorComponentValue(q, j, l);
A_c(i, j, k, l)
+= w * eProblem.TractionStiffness(Time(), faceElem.QPoint(q),E.Bnd(f))
* Nphi_ik * phi_jl;
}
}
}
}
}
}
}
}
}
}
void DGElasticity::Project(const Vector &fine, Vector &u) const { void DGElasticity::Project(const Vector &fine, Vector &u) const {
int dlevel = fine.SpaceLevel() - u.SpaceLevel(); int dlevel = fine.SpaceLevel() - u.SpaceLevel();
...@@ -854,14 +892,11 @@ double DGElasticity::H1Error(const Vector &u, const Vector &reference) const { ...@@ -854,14 +892,11 @@ double DGElasticity::H1Error(const Vector &u, const Vector &reference) const {
} }
void DGElasticity::prestress(const Vector &u, Vector &pStress) { void DGElasticity::prestress(const Vector &u, Vector &pStress) {
Prestress = std::make_unique<Vector>(0.0, u);
//TODO THis has to be updated
for (cell c = u.cells(); c != u.cells_end(); ++c) { for (cell c = u.cells(); c != u.cells_end(); ++c) {
const auto &cellMat = eProblem.GetMaterial(c); const auto &cellMat = eProblem.GetMaterial(c);
DGVectorFieldElement E(u, *c); DGVectorFieldElement E(u, *c);
DGSizedRowBndValues r_c(*Prestress, *c, E.shape_size()); DGSizedRowBndValues r_c(pStress, *c, E.shape_size());
Tensor Q = OrientationMatrix(c.GetData()); Tensor Q = OrientationMatrix(c.GetData());
...@@ -880,10 +915,119 @@ void DGElasticity::prestress(const Vector &u, Vector &pStress) { ...@@ -880,10 +915,119 @@ void DGElasticity::prestress(const Vector &u, Vector &pStress) {
} }
} }
} }
} for (int f = 0; f < c.Faces(); ++f) {
double s = cellMat.Isochoric().GetMaxParameter() *
(pow(degree, 2) * penalty) / u.GetMesh().MaxMeshWidth();
DGVectorFieldFaceElement FE(u, *c, f);
if (E.Bnd(f) == 1) { // Dirichlet-RB
for (int q = 0; q < FE.nQ(); ++q) {
double w = FE.QWeight(q);
Point z = FE.QPoint(q);
VectorField N = FE.QNormal(q);
VectorField U = eProblem.Deformation(Time(), z);
VectorField u_diff = FE.VectorValue(q, u) - U;
Tensor F = DeformationGradient(FE, q, u);
auto Fiso = eProblem.IsochoricPart(F);
for (int i = 0; i < FE.shape_size(); ++i)
for (int k = 0; k < E.Dim(); ++k) {
auto u_ik = FE.VectorComponentValue(q, i, k);
auto Du_ik = FE.VectorRowGradient(q, i, k);
double SN = cellMat.TotalDerivative(Fiso, F, Q, Product(u_ik, N));
double SN_ik = cellMat.TotalSecondDerivative(Fiso, F, Q, Du_ik, Product(u_diff, N));
r_c(i, k) -= w * (SN + SN_ik * sign - s * u_diff * u_ik);
}
}
} else if (E.Bnd(f) == 199) { // Dirichlet-RB
for (int q = 0; q < FE.nQ(); ++q) {
double w = FE.QWeight(q);
Point z = FE.QPoint(q);
VectorField N = FE.QNormal(q);
VectorField u_diff = FE.VectorValue(q, u);
Tensor F = DeformationGradient(FE, q, u);
auto Fiso = eProblem.IsochoricPart(F);
for (int i = 0; i < FE.shape_size(); ++i)
for (int k = 0; k < E.Dim(); ++k) {
auto u_ik = FE.VectorComponentValue(q, i, k);
auto Du_ik = FE.VectorRowGradient(q, i, k);
double SN = cellMat.TotalDerivative(Fiso, F, Q, Product(u_ik, N));
double SN_ik = cellMat.TotalSecondDerivative(Fiso, F, Q, Du_ik, Product(u_diff, N));
r_c(i, k) -= w * (SN + SN_ik * sign - s * u_diff * u_ik);
}
}
} else if (E.Bnd(f) == -1) {
cell cf = pStress.find_neighbour_cell(c, f);
if (cf() < c()) continue;
int f1 = pStress.find_neighbour_face_id(c.Face(f), cf);
DGVectorFieldElement E_nghbr(u, *cf);
DGVectorFieldFaceElement FE_nghbr(pStress, *cf, f1);
DGSizedRowBndValues r_cf(pStress, *cf, E_nghbr.shape_size());
s = cellMat.Isochoric().GetMaxParameter() * (pow(degree, 2) * penalty) /
u.GetMesh().MaxMeshWidth();
for (int q = 0; q < FE.nQ(); ++q) {
double w = FE.QWeight(q);
Point z = FE.QPoint(q);
VectorField N = FE.QNormal(q);
// für erstes Face-Element
VectorField U = FE.VectorValue(q, u);
Tensor F = DeformationGradient(FE, q, u);
auto Fiso = eProblem.IsochoricPart(F);
Tensor inverseFA = mat::active_deformation(Q, FE.VectorValue(q, *gamma));
F = F * inverseFA;
// das gleiche für das zweite Face-Element
int q1 = FE.findQPointID(FE_nghbr, z);
VectorField u_nghbr = FE_nghbr.VectorValue(q1, u);
Tensor F_nghbr = DeformationGradient(FE_nghbr, q1, u);
auto Fiso_nghbr = eProblem.IsochoricPart(F_nghbr);
Tensor inverseFA_nghbr = mat::active_deformation(Q, FE_nghbr.VectorValue(q1, *gamma));
F_nghbr = F_nghbr * inverseFA_nghbr;
for (int i = 0; i < FE.shape_size(); ++i) {
for (int k = 0; k < E.Dim(); ++k) {
auto u_ik = FE.VectorComponentValue(q, i, k);
auto Du_ik = FE.VectorRowGradient(q, i, k);
Prestress->ClearDirichletValues(); double SN = cellMat.TotalDerivative(Fiso, F, Q, Product(u_ik, N));
Prestress->Collect(); double SN_ik = cellMat.TotalSecondDerivative(Fiso, F, Q, Du_ik, Product(U, N));
double SN_nghbr = cellMat.TotalDerivative(Fiso_nghbr, F_nghbr, Q, Product(u_ik, N));
double SN_ik_nghbr = cellMat.TotalSecondDerivative(Fiso, F, Q, Du_ik,
Product(u_nghbr, N));
r_c(i, k) += w * (s * U * u_ik - 0.5 * SN - 0.5 * SN_ik * sign);
r_c(i, k) += w * (-s * u_nghbr * u_ik - 0.5 * SN_nghbr + 0.5 * SN_ik_nghbr * sign);
}
}
for (int j = 0; j < FE_nghbr.shape_size(); ++j)
for (int k = 0; k < E.Dim(); ++k) {
auto u_jk = FE_nghbr.VectorComponentValue(q1, j, k);
auto Du_jk = FE_nghbr.VectorRowGradient(q1, j, k);
double SN = cellMat.TotalDerivative(Fiso, F, Q, Product(u_jk, N));
double SN_jk = cellMat.TotalSecondDerivative(Fiso_nghbr, F_nghbr, Q, Du_jk,
Product(U, N));
double SN_nghbr = cellMat.TotalDerivative(Fiso_nghbr, F_nghbr, Q, Product(u_jk, N));
double SN_jk_nghbr = cellMat.TotalSecondDerivative(Fiso_nghbr, F_nghbr, Q, Du_jk,
Product(u_nghbr, N));
r_cf(j, k) += w * (s * u_nghbr * u_jk + 0.5 * SN_nghbr + 0.5 * SN_jk_nghbr * sign);
r_cf(j, k) += w * (-s * U * u_jk + 0.5 * SN - 0.5 * SN_jk * sign);
}
}
}
}
}
} }
...@@ -1278,3 +1422,8 @@ std::pair<double, double> DGElasticity::detF(const Vector &u) const { ...@@ -1278,3 +1422,8 @@ std::pair<double, double> DGElasticity::detF(const Vector &u) const {
return {Jmin, Jmax}; return {Jmin, Jmax};
} }
void DGElasticity::SetInitialValue(Vector &u) {
IElasticity::SetInitialValue(u);
}
\ No newline at end of file
src/elasticity/assemble/DGElasticity.hpp
View file @ 8f9a221d
...@@ -39,6 +39,8 @@ public: ...@@ -39,6 +39,8 @@ public:
void Jacobi(const cell &c, const Vector &u, Matrix &A) const override; void Jacobi(const cell &c, const Vector &u, Matrix &A) const override;
void TractionStiffness(Matrix &matrix) const override;
void MassMatrix(Matrix &massMatrix) const override; void MassMatrix(Matrix &massMatrix) const override;
void SystemMatrix(Matrix &systemMatrix) const override; void SystemMatrix(Matrix &systemMatrix) const override;
...@@ -129,6 +131,8 @@ public: ...@@ -129,6 +131,8 @@ public:
} }
std::pair<double, double> detF(const Vector &u) const override; std::pair<double, double> detF(const Vector &u) const override;
void SetInitialValue(Vector &u) override;
}; };
#endif //DGELASTICITY_HPP #endif //DGELASTICITY_HPP
src/elasticity/assemble/EGElasticity.cpp
View file @ 8f9a221d
...@@ -344,7 +344,7 @@ void EGElasticity::Jacobi(const cell &c, const Vector &U, Matrix &A) const { ...@@ -344,7 +344,7 @@ void EGElasticity::Jacobi(const cell &c, const Vector &U, Matrix &A) const {
auto Fiso = eProblem.IsochoricPart(F); auto Fiso = eProblem.IsochoricPart(F);
Tensor inverseFA = mat::active_deformation(Q, E.VectorValue(q, *gamma)); Tensor inverseFA = mat::active_deformation(Q, E.VectorValue(q, *gamma));
F = F * inverseFA; F = F * inverseFA;
double T = eProblem.ActiveStress(Time(), E.QPoint(q)); // double T = eProblem.ActiveStress(Time(), E.QPoint(q));
auto theta = E.PenaltyVectorValue(q); auto theta = E.PenaltyVectorValue(q);
auto Dtheta = E.PenaltyVectorGradient(q) * transpose(inverseFA); auto Dtheta = E.PenaltyVectorGradient(q) * transpose(inverseFA);
...@@ -586,6 +586,56 @@ void EGElasticity::Jacobi(const cell &c, const Vector &U, Matrix &A) const { ...@@ -586,6 +586,56 @@ void EGElasticity::Jacobi(const cell &c, const Vector &U, Matrix &A) const {
} }
} }
void EGElasticity::TractionStiffness(Matrix &matrix) const {
matrix = 0;
for (cell c = matrix.cells(); c != matrix.cells_end(); ++c) {
MixedEGVectorFieldElement E(matrix, *c);
MixedRowEntries A_c(matrix, *c);
for (int f = 0; f < c.Faces(); ++f) {
MixedEGVectorFieldFaceElement faceElem(matrix, *c, f);
if (matrix.OnBoundary(*c, f)) {
int bc = matrix.GetMesh().BoundaryFacePart(c.Face(f));
if ((bc >= ROBIN_BC_EPI) && (bc <= ROBIN_BC_BASE)) {
for (int q = 0; q < faceElem.nQ(); ++q) {
VectorField n = faceElem.QNormal(q);
auto N = Product(n, n);
double w = faceElem.QWeight(q);
auto theta = faceElem.PenaltyVectorValue(q);
auto Ntheta = N * faceElem.PenaltyVectorValue(q);
for (int i = 0; i < faceElem.ValueSize(); ++i) {
for (int k = 0; k < faceElem.Dim(); ++k) {
auto phi_ik = faceElem.VectorComponentValue(q, i, k);
auto Nphi_ik = N * faceElem.VectorComponentValue(q, i, k);
for (int j = 0; j < faceElem.ValueSize(); ++j) {
for (int l = 0; l < faceElem.Dim(); ++l) {
auto phi_jl = faceElem.VectorComponentValue(q, j, l);
A_c(faceElem.ValueIndex(), faceElem.ValueIndex(),i, j, k, l)
+= w * eProblem.TractionStiffness(Time(), faceElem.QPoint(q),E.Bnd(f))
* Nphi_ik * phi_jl;
}
}
A_c(faceElem.ValueIndex(), faceElem.PenaltyIndex(),i, 0, k, 0)
+= w * eProblem.TractionStiffness(Time(), faceElem.QPoint(q),E.Bnd(f))
* Ntheta * phi_ik;
A_c(faceElem.PenaltyIndex(), faceElem.ValueIndex(),0, i, 0, k)
+= w * eProblem.TractionStiffness(Time(), faceElem.QPoint(q),E.Bnd(f))
* Nphi_ik * theta;
}
}
A_c(faceElem.PenaltyIndex(), faceElem.PenaltyIndex(),0, 0, 0, 0)
+= w * eProblem.TractionStiffness(Time(), faceElem.QPoint(q),E.Bnd(f))
* Ntheta * theta;
}
}
}
}
}
}
double EGElasticity::L2Norm(const Vector &u) const { double EGElasticity::L2Norm(const Vector &u) const {
double err = 0; double err = 0;
for (cell c = u.cells(); c != u.cells_end(); ++c) { for (cell c = u.cells(); c != u.cells_end(); ++c) {
...@@ -743,14 +793,14 @@ double EGElasticity::EnergyError(const Vector &u) const { ...@@ -743,14 +793,14 @@ double EGElasticity::EnergyError(const Vector &u) const {
void EGElasticity::prestress(const Vector &u, Vector &pStress) { void EGElasticity::prestress(const Vector &u, Vector &pStress) {
Prestress = std::make_unique<Vector>(0.0, u);
for (cell c = u.cells(); c != u.cells_end(); ++c) { for (cell c = u.cells(); c != u.cells_end(); ++c) {
const auto &cellMat = eProblem.GetMaterial(c); const auto &cellMat = eProblem.GetMaterial(c);
MixedEGVectorFieldElement E(u, *c); MixedEGVectorFieldElement E(u, *c);
MixedRowValues r_c(*Prestress, *c); MixedRowValues r_c(pStress, *c);
Tensor Q = OrientationMatrix(c.GetData()); Tensor Q = OrientationMatrix(c.GetData());
for (int q = 0; q < E.nQ(); ++q) { for (int q = 0; q < E.nQ(); ++q) {
double w = E.QWeight(q); double w = E.QWeight(q);
...@@ -759,7 +809,7 @@ void EGElasticity::prestress(const Vector &u, Vector &pStress) { ...@@ -759,7 +809,7 @@ void EGElasticity::prestress(const Vector &u, Vector &pStress) {
for (int i = 0; i < E.ValueSize(); ++i) { for (int i = 0; i < E.ValueSize(); ++i) {
for (int k = 0; k < E.Dim(); ++k) { for (int k = 0; k < E.Dim(); ++k) {
auto F_ik = Q * E.VectorRowGradient(q, i, k); auto F_ik = E.VectorRowGradient(q, i, k);
// Passive material response // Passive material response
r_c(E.ValueIndex(), i, k) += w * cellMat.TotalDerivative(Fiso, F, Q, F_ik); r_c(E.ValueIndex(), i, k) += w * cellMat.TotalDerivative(Fiso, F, Q, F_ik);
...@@ -768,9 +818,153 @@ void EGElasticity::prestress(const Vector &u, Vector &pStress) { ...@@ -768,9 +818,153 @@ void EGElasticity::prestress(const Vector &u, Vector &pStress) {
auto Dtheta = E.PenaltyVectorGradient(q); auto Dtheta = E.PenaltyVectorGradient(q);
r_c(E.PenaltyIndex(), 0, 0) += w * cellMat.TotalDerivative(Fiso, F, Q, Dtheta); r_c(E.PenaltyIndex(), 0, 0) += w * cellMat.TotalDerivative(Fiso, F, Q, Dtheta);
} }
for (int f = 0; f < c.Faces(); ++f) {
double s = cellMat.Isochoric().GetMaxParameter() *
(pow(degree, 2) * penalty) / u.GetMesh().MaxMeshWidth();
MixedEGVectorFieldFaceElement FE(u, *c, f);
if (E.Bnd(f) == 1) {
for (int q = 0; q < FE.nQ(); ++q) {
double w = FE.QWeight(q);
VectorField N = FE.QNormal(q);
VectorField u_diff = FE.VectorValue(q, u)- eProblem.Deformation(Time(), FE.QPoint(q));
Tensor F = DeformationGradient(FE, q, u);
auto Fiso = eProblem.IsochoricPart(F);
Tensor inverseFA = mat::active_deformation(Q, FE.VectorValue(q, *gamma));
F = F * inverseFA;
for (int i = 0; i < FE.ValueSize(); ++i) {
for (int k = 0; k < FE.Dim(); ++k) {
auto u_ik = FE.VectorComponentValue(q, i, k);
auto Du_ik = FE.VectorRowGradient(q, i, k);
double SN = cellMat.TotalDerivative(Fiso, F, Q, Product(u_ik, N));
double SN_ik = cellMat.TotalSecondDerivative(Fiso, F, Q, Du_ik,
Product(u_diff, N));
r_c(FE.ValueIndex(), i, k) -= w * (SN_ik * sign + SN - s * u_diff * u_ik);
}
}
auto theta = FE.PenaltyVectorValue(q);
auto Dtheta = FE.PenaltyVectorGradient(q);
double SN_theta = cellMat.TotalDerivative(Fiso, F, Q, Product(theta, N));
double SN_ik_Dtheta = cellMat.TotalSecondDerivative(Fiso, F, Q, Dtheta,
Product(u_diff, N));
r_c(FE.PenaltyIndex(), 0, 0) -= w * (SN_ik_Dtheta * sign + SN_theta - s * u_diff * theta);
}
} else if( E.Bnd(f) == 199) {
for (int q = 0; q < FE.nQ(); ++q) {
double w = FE.QWeight(q);
VectorField N = FE.QNormal(q);
VectorField u_diff = FE.VectorValue(q, u);
Tensor F = DeformationGradient(FE, q, u);
auto Fiso = eProblem.IsochoricPart(F);
Tensor inverseFA = mat::active_deformation(Q, FE.VectorValue(q, *gamma));
F = F * inverseFA;
for (int i = 0; i < FE.ValueSize(); ++i) {
for (int k = 0; k < FE.Dim(); ++k) {
auto u_ik = FE.VectorComponentValue(q, i, k);
auto Du_ik = FE.VectorRowGradient(q, i, k);
double SN = cellMat.TotalDerivative(Fiso, F, Q, Product(u_ik, N));
double SN_ik = cellMat.TotalSecondDerivative(Fiso, F, Q, Du_ik,
Product(u_diff, N));
r_c(FE.ValueIndex(), i, k) -= w * ( SN_ik * sign + SN - s * u_diff * u_ik);
}
}
auto theta = FE.PenaltyVectorValue(q);
auto Dtheta = FE.PenaltyVectorGradient(q);
double SN_theta = cellMat.TotalDerivative(Fiso, F, Q, Product(theta, N));
double SN_Dtheta = cellMat.TotalSecondDerivative(Fiso, F, Q, Dtheta,
Product(u_diff, N));
r_c(FE.PenaltyIndex(), 0, 0) -= w *(SN_Dtheta * sign + SN_theta - s * u_diff * theta);
}
} else if (E.Bnd(f) == -1) {
cell cf = pStress.find_neighbour_cell(c, f);
if (cf() < c()) continue;
int f1 = pStress.find_neighbour_face_id(c.Face(f), cf);
MixedEGVectorFieldElement E_nghbr(u, *cf);
MixedEGVectorFieldFaceElement FE_nghbr(pStress, *cf, f1);
MixedRowValues r_cf(pStress, *cf);
for (int q = 0; q < FE.nQ(); ++q) {
double w = FE.QWeight(q);
Point z = FE.QPoint(q);
VectorField N = FE.QNormal(q);
VectorField U = FE.VectorValue(q, u);
Tensor F = DeformationGradient(FE, q, u);
auto Fiso = eProblem.IsochoricPart(F);
int q1 = FE.FindQPointID(FE_nghbr, z);
VectorField u_nghbr = FE_nghbr.VectorValue(q1, u);
Tensor F_nghbr = DeformationGradient(FE_nghbr, q1, u);
auto Fiso_nghbr = eProblem.IsochoricPart(F_nghbr);
for (int i = 0; i < FE.ValueSize(); ++i) {
for (int k = 0; k < FE.Dim(); ++k) {
auto u_ik = FE.VectorComponentValue(q, i, k);
auto Du_ik = FE.VectorRowGradient(q, i, k);
double SN = cellMat.TotalDerivative(Fiso, F, Q, Product(u_ik, N));
double SN_ik = cellMat.TotalSecondDerivative(Fiso, F, Q, Du_ik, Product(U, N));
double SN_nghbr = cellMat.TotalDerivative(Fiso_nghbr, F_nghbr, Q, Product(u_ik, N));
double SN_ik_nghbr = cellMat.TotalSecondDerivative(Fiso, F, Q, Du_ik,
Product(u_nghbr, N));
r_c(FE.ValueIndex(), i, k) += w * (s * U * u_ik - 0.5 * SN - 0.5 * SN_ik * sign);
r_c(FE.ValueIndex(), i, k) +=
w * (-s * u_nghbr * u_ik - 0.5 * SN_nghbr + 0.5 * SN_ik_nghbr * sign);
}
}
for (int j = 0; j < FE_nghbr.ValueSize(); ++j) {
for (int l = 0; l < FE_nghbr.Dim(); ++l) {
auto u_jl = FE_nghbr.VectorComponentValue(q1, j, l);
auto Du_jl = FE_nghbr.VectorRowGradient(q1, j, l);
double SN = cellMat.TotalDerivative(Fiso, F, Q, Product(u_jl, N));
double SN_jl = cellMat.TotalSecondDerivative(Fiso_nghbr, F_nghbr, Q, Du_jl,
Product(U, N));
double SN_nghbr = cellMat.TotalDerivative(Fiso_nghbr, F_nghbr, Q, Product(u_jl, N));
double SN_jl_nghbr = cellMat.TotalSecondDerivative(Fiso_nghbr, F_nghbr, Q, Du_jl,
Product(u_nghbr, N));
r_cf(FE_nghbr.ValueIndex(), j, l) +=
w * (s * u_nghbr * u_jl + 0.5 * SN_nghbr + 0.5 * SN_jl_nghbr * sign);
r_cf(FE_nghbr.ValueIndex(), j, l) +=
w * (-s * U * u_jl + 0.5 * SN - 0.5 * SN_jl * sign);
}
}
auto theta = FE.PenaltyVectorValue(q);
auto Dtheta = FE.PenaltyVectorGradient(q);
double SN_theta = cellMat.TotalDerivative(Fiso, F, Q, Product(theta, N));
double SN_ik_Dtheta = cellMat.TotalSecondDerivative(Fiso, F, Q, Dtheta, Product(U, N));
double SN_nghbr_theta = cellMat.TotalDerivative(Fiso_nghbr, F_nghbr, Q,
Product(theta, N));
double SN_ik_nghbr_Dtheta = cellMat.TotalSecondDerivative(Fiso, F, Q, Dtheta,
Product(u_nghbr, N));
r_c(FE.PenaltyIndex(), 0, 0) +=
w * (s * U * theta - 0.5 * SN_theta - 0.5 * SN_ik_Dtheta * sign);
r_c(FE.PenaltyIndex(), 0, 0) +=
w * (-s * u_nghbr * theta - 0.5 * SN_nghbr_theta + 0.5 * SN_ik_nghbr_Dtheta * sign);
auto theta_nghbr = FE_nghbr.PenaltyVectorValue(q1);
auto Dtheta_nghbr = FE_nghbr.PenaltyVectorGradient(q1);
double SN_theta_nghbr = cellMat.TotalDerivative(Fiso, F, Q, Product(theta_nghbr, N));
double SN_jl_Dtheta_nghbr = cellMat.TotalSecondDerivative(Fiso_nghbr, F_nghbr, Q,
Dtheta_nghbr,
Product(U, N));
double SN_nghbr_theta_nghbr = cellMat.TotalDerivative(Fiso_nghbr, F_nghbr, Q,
Product(theta_nghbr, N));
double SN_jl_nghbr_Dtheta_nghbr = cellMat.TotalSecondDerivative(Fiso_nghbr, F_nghbr, Q,
Dtheta_nghbr,
Product(u_nghbr, N));
r_cf(FE_nghbr.PenaltyIndex(), 0, 0) += w * (s * u_nghbr * theta_nghbr +
0.5 * SN_nghbr_theta_nghbr +
0.5 * SN_jl_nghbr_Dtheta_nghbr * sign);//?
r_cf(FE_nghbr.PenaltyIndex(), 0, 0) +=
w * (-s * U * theta_nghbr + 0.5 * SN_theta_nghbr - 0.5 * SN_jl_Dtheta_nghbr * sign);
}
}
}
} }
Prestress->ClearDirichletValues();
Prestress->Collect();
} }
double EGElasticity::PressureError(const Vector &u) const { double EGElasticity::PressureError(const Vector &u) const {
...@@ -1244,3 +1438,6 @@ std::pair<double, double> EGElasticity::detF(const Vector &u) const { ...@@ -1244,3 +1438,6 @@ std::pair<double, double> EGElasticity::detF(const Vector &u) const {
return {Jmin, Jmax}; return {Jmin, Jmax};
} }
void EGElasticity::SetInitialValue(Vector &u) {
IElasticity::SetInitialValue(u);
}
src/elasticity/assemble/EGElasticity.hpp
View file @ 8f9a221d
...@@ -54,6 +54,8 @@ public: ...@@ -54,6 +54,8 @@ public:
void Jacobi(const cell &c, const Vector &u, Matrix &A) const override; void Jacobi(const cell &c, const Vector &u, Matrix &A) const override;
void TractionStiffness(Matrix &matrix) const override;
void MassMatrix(Matrix &massMatrix) const override; void MassMatrix(Matrix &massMatrix) const override;
void SystemMatrix(Matrix &systemMatrix) const override; void SystemMatrix(Matrix &systemMatrix) const override;
...@@ -112,7 +114,14 @@ public: ...@@ -112,7 +114,14 @@ public:
std::array<double, 4> ChamberVolume(const Vector &u) const override; std::array<double, 4> ChamberVolume(const Vector &u) const override;
std::pair<double, double> detF(const Vector &u) const override; std::pair<double, double> detF(const Vector &u) const override;
void SetInitialValue(Vector &u) override;
void Scale (double s) const {
Material &cellMat = eProblem.GetMaterial();
//cellMat.ScalePenalty(s*s);
}
}; };
......
src/elasticity/assemble/IElasticity.hpp
View file @ 8f9a221d
...@@ -81,8 +81,8 @@ public: ...@@ -81,8 +81,8 @@ public:
using IAssemble::Residual; using IAssemble::Residual;
double Residual(const Vector &u, Vector &r) const override { double Residual(const Vector &u, Vector &r) const override {
r=0;
r = *Prestress; r = *Prestress;
r = 0;
for (cell ce = u.cells(); ce != u.cells_end(); ++ce) { for (cell ce = u.cells(); ce != u.cells_end(); ++ce) {
Residual(ce, u, r); Residual(ce, u, r);
} }
......
src/elasticity/solvers/IterativePressureSolver.cpp
View file @ 8f9a221d
...@@ -21,7 +21,6 @@ bool IterativePressureSolver::Method(IElasticity &assemble, Vector &u) { ...@@ -21,7 +21,6 @@ bool IterativePressureSolver::Method(IElasticity &assemble, Vector &u) {
<< endl; << endl;
Vector uNew(u); Vector uNew(u);
Initialize(assemble, uNew); Initialize(assemble, uNew);
while (!iteration.IsFinished()) { while (!iteration.IsFinished()) {
bool converged = Step(assemble, uNew); bool converged = Step(assemble, uNew);
...@@ -52,7 +51,7 @@ bool IterativePressureSolver::Step(IElasticity &assemble, Vector &u) { ...@@ -52,7 +51,7 @@ bool IterativePressureSolver::Step(IElasticity &assemble, Vector &u) {
Config::Get("WithPrestress", prestress); Config::Get("WithPrestress", prestress);
std::string activeDeformation; std::string activeDeformation;
Config::Get("activeDeformation", activeDeformation); Config::Get("activeDeformation", activeDeformation);
if (activeDeformation != "" && prestress == false) { if (activeDeformation != "" && !prestress) {
assemble.UpdateStretch(u); assemble.UpdateStretch(u);
} }
iteration.NextTimeStep(); iteration.NextTimeStep();
...@@ -84,17 +83,29 @@ bool KlotzPressureSolver::Step(IElasticity &assemble, Vector &u) { ...@@ -84,17 +83,29 @@ bool KlotzPressureSolver::Step(IElasticity &assemble, Vector &u) {
return conv; return conv;
} }
//TODO implement EG Prolongarte and use this newer version again!!!
//void CalculatePrestress(IElasticity &assemble, Vector &u) {
// int prestressSteps{1};
// Config::Get("PrestressSteps", prestressSteps);
// int prestressLevel=0;
// Config::Get("PrestressLevel", prestressLevel);
// Vector uPrestressLevel(assemble.GetSharedDisc(),prestressLevel);
// uPrestressLevel.Accumulate();
// IterativePressureSolver prestressSolver(assemble.GetElasticityProblem(), prestressSteps);
// prestressSolver.Method(assemble, uPrestressLevel);
// LagrangeTransfer lagrangeT(uPrestressLevel,u);
// lagrangeT.Prolongate(uPrestressLevel,u);
// assemble.InitializePrestress(u);
// //u.Clear();
//}
//old version
void CalculatePrestress(IElasticity &assemble, Vector &u) { void CalculatePrestress(IElasticity &assemble, Vector &u) {
int prestressSteps{1}; int prestressSteps{1};
Config::Get("PrestressSteps", prestressSteps); Config::Get("PrestressSteps", prestressSteps);
int prestressLevel=0;
Config::Get("PrestressLevel", prestressLevel);
Vector uPrestressLevel(assemble.GetSharedDisc(),prestressLevel);
uPrestressLevel.Accumulate();
IterativePressureSolver prestressSolver(assemble.GetElasticityProblem(), prestressSteps); IterativePressureSolver prestressSolver(assemble.GetElasticityProblem(), prestressSteps);
prestressSolver.Method(assemble, uPrestressLevel); prestressSolver.Method(assemble, u);
LagrangeTransfer lagrangeT(uPrestressLevel,u);
lagrangeT.Prolongate(uPrestressLevel,u);
assemble.InitializePrestress(u); assemble.InitializePrestress(u);
//u.Clear(); //u.Clear();
} }
\ No newline at end of file
src/elasticity/solvers/StaticSolver.cpp
View file @ 8f9a221d
...@@ -9,6 +9,8 @@ void StaticSolver::Initialize(const IElasticity &assemble, Vector &u) { ...@@ -9,6 +9,8 @@ void StaticSolver::Initialize(const IElasticity &assemble, Vector &u) {
double StaticSolver::calculateResidualUpdate(const IElasticity &assemble, const Vector &u, double StaticSolver::calculateResidualUpdate(const IElasticity &assemble, const Vector &u,
Vector &defect) const { Vector &defect) const {
defect = 0;
assemble.Residual(u, defect); assemble.Residual(u, defect);
defect += *residualMatrix * u; defect += *residualMatrix * u;
......
test/elasticity/CMakeLists.txt
View file @ 8f9a221d
...@@ -19,7 +19,7 @@ add_mpp_test(TestLaplaceElasticity ELASTICITY) ...@@ -19,7 +19,7 @@ add_mpp_test(TestLaplaceElasticity ELASTICITY)
#add_mpp_test(TestDynamicBoundary ELASTICITY) #add_mpp_test(TestDynamicBoundary ELASTICITY)
#add_mpp_test(TestOscillation ELASTICITY) #add_mpp_test(TestOscillation ELASTICITY)
#add_mpp_test(TestCFLCondition ELASTICITY) #add_mpp_test(TestCFLCondition ELASTICITY)
#add_mpp_test(TestPrestress ELASTICITY) add_mpp_test(TestPrestress ELASTICITY)
add_mpp_test(TestVolume ELASTICITY) add_mpp_test(TestVolume ELASTICITY)
#add_mpp_test(TestVolumePenalty ELASTICITY) #add_mpp_test(TestVolumePenalty ELASTICITY)
...@@ -36,7 +36,7 @@ add_mpi_test(TestQuadraticBeamProblem ELASTICITY) ...@@ -36,7 +36,7 @@ add_mpi_test(TestQuadraticBeamProblem ELASTICITY)
#add_mpi_test(TestElasticityBlock ELASTICITY) #add_mpi_test(TestElasticityBlock ELASTICITY)
add_mpi_test(TestDynamicBoundary ELASTICITY) add_mpi_test(TestDynamicBoundary ELASTICITY)
#add_mpi_test(TestCFLCondition ELASTICITY) #add_mpi_test(TestCFLCondition ELASTICITY)
#add_mpi_test(TestPrestress ELASTICITY) add_mpi_test(TestPrestress ELASTICITY)
add_mpi_test(TestVolume ELASTICITY) add_mpi_test(TestVolume ELASTICITY)
#add_mpi_test(TestOscillation ELASTICITY) #add_mpi_test(TestOscillation ELASTICITY)
test/elasticity/TestPrestress.cpp
View file @ 8f9a221d
...@@ -27,13 +27,14 @@ protected: ...@@ -27,13 +27,14 @@ protected:
testConfig["MechDynamics"] = "Static"; testConfig["MechDynamics"] = "Static";
testConfig["MechDiscretization"] = GetParam().discType; testConfig["MechDiscretization"] = GetParam().discType;
testConfig["PressureSteps"] = "1"; testConfig["PressureSteps"] = "1";
testConfig["DGSign"] = "-1";
testConfig["DGPenalty"] = "10";
testConfig["PrestressSteps"] = "5"; testConfig["PrestressSteps"] = "5";
testConfig["WithPrestress"] = "true"; testConfig["WithPrestress"] = "true";
testConfig["MechEpsilon"] = "1e-12"; testConfig["MechEpsilon"] = "1e-12";
testConfig["NewtonEpsilon"] = "1e-11"; testConfig["NewtonEpsilon"] = "1e-10";
testConfig["NewtonReduction"] = "1e-12"; testConfig["NewtonReduction"] = "1e-11";
testConfig["WithPrestress"] = "true";
//testConfig["Mesh"] = GetParam().meshName; //testConfig["Mesh"] = GetParam().meshName;
...@@ -69,10 +70,14 @@ TEST_P(PrestressTest, TestInitialValue) { ...@@ -69,10 +70,14 @@ TEST_P(PrestressTest, TestInitialValue) {
mout << "solution " << solution << endl; mout << "solution " << solution << endl;
mout << "num " << cmMain->EvaluateQuantity(solution, "L2") << endl << endl; mout << "num " << cmMain->EvaluateQuantity(solution, "L2") << endl << endl;
for (row r = solution.rows(); r != solution.rows_end(); ++r) { for (row r = solution.rows(); r != solution.rows_end(); ++r) {
for (int i = 0; i < r.n(); ++i) { for (int i = 0; i < r.NumberOfDofs(); ++i) {
EXPECT_NEAR(solution(r, i), 0.0, NEWTON_TOLERANCE); EXPECT_NEAR(solution(r, i), 0.0, NEWTON_TOLERANCE);
} }
} }
Vector exactSolution(cmMain->GetDisplacement());
cmMain->SetDisplacement(true, exactSolution);
mout << "exact " << exactSolution << endl << endl;
} }
...@@ -84,7 +89,7 @@ INSTANTIATE_TEST_SUITE_P(ConformingBeam, PrestressTest, Values( ...@@ -84,7 +89,7 @@ INSTANTIATE_TEST_SUITE_P(ConformingBeam, PrestressTest, Values(
TestPrestressParameter({0, 2, {10.0, 0.0, 1.0, 0.0}, "PrestressBeam", "Conforming"}), TestPrestressParameter({0, 2, {10.0, 0.0, 1.0, 0.0}, "PrestressBeam", "Conforming"}),
TestPrestressParameter({1, 2, {10.0, 0.0, 1.0, 0.0}, "PrestressBeam", "Conforming"}) TestPrestressParameter({1, 2, {10.0, 0.0, 1.0, 0.0}, "PrestressBeam", "Conforming"})
)); ));
/*
INSTANTIATE_TEST_SUITE_P(DGBeam, PrestressTest, Values( INSTANTIATE_TEST_SUITE_P(DGBeam, PrestressTest, Values(
//TestPrestressParameter({0, 0, {0.0, 0.0,-0.0005,0.0}, "BenchmarkBeam3DTet"}), //TestPrestressParameter({0, 0, {0.0, 0.0,-0.0005,0.0}, "BenchmarkBeam3DTet"}),
TestPrestressParameter({0, 1, {10.0, 0.0, 1.0, 0.0}, "PrestressBeam", "DG"}), TestPrestressParameter({0, 1, {10.0, 0.0, 1.0, 0.0}, "PrestressBeam", "DG"}),
...@@ -102,7 +107,7 @@ INSTANTIATE_TEST_SUITE_P(EGBeam, PrestressTest, Values( ...@@ -102,7 +107,7 @@ INSTANTIATE_TEST_SUITE_P(EGBeam, PrestressTest, Values(
TestPrestressParameter({0, 2, {10.0, 0.0, 1.0, 0.0}, "PrestressBeam", "EG"}), TestPrestressParameter({0, 2, {10.0, 0.0, 1.0, 0.0}, "PrestressBeam", "EG"}),
TestPrestressParameter({1, 2, {10.0, 0.0, 1.0, 0.0}, "PrestressBeam", "EG"}) TestPrestressParameter({1, 2, {10.0, 0.0, 1.0, 0.0}, "PrestressBeam", "EG"})
)); ));
*/
//TODO Why won't this work? //TODO Why won't this work?
/* /*
...@@ -117,4 +122,4 @@ INSTANTIATE_TEST_SUITE_P(OnEllipsoid, PrestressTest, Values( ...@@ -117,4 +122,4 @@ INSTANTIATE_TEST_SUITE_P(OnEllipsoid, PrestressTest, Values(
int main(int argc, char **argv) { int main(int argc, char **argv) {
MppTest mppTest = MppTestBuilder(argc, argv).WithPPM().WithoutDefaultConfig().WithScreenLogging().WithFileLogging(); MppTest mppTest = MppTestBuilder(argc, argv).WithPPM().WithoutDefaultConfig().WithScreenLogging().WithFileLogging();
return mppTest.RUN_ALL_MPP_TESTS(); return mppTest.RUN_ALL_MPP_TESTS();
} }
\ No newline at end of file

Help | Imprint | Privacy policy | Accessibility | Contact