From 154533e332c02da7087978bd641377aed4553cd6 Mon Sep 17 00:00:00 2001
From: "jonathan.froehlich" <jonathan.froehlich@kit.edu>
Date: Mon, 28 Feb 2022 19:08:41 +0100
Subject: [PATCH] Fixed DG Projection
---
.../src/elasticity/assemble/DGElasticity.cpp | 32 +++++++++++++------
.../problem/ElasticityBenchmarkProblem.hpp | 2 +-
tools/jobsystem/run_dg_benchmarks.py | 0
3 files changed, 24 insertions(+), 10 deletions(-)
create mode 100644 tools/jobsystem/run_dg_benchmarks.py
diff --git a/cardmech/src/elasticity/assemble/DGElasticity.cpp b/cardmech/src/elasticity/assemble/DGElasticity.cpp
index b49274155..526b16020 100644
--- a/cardmech/src/elasticity/assemble/DGElasticity.cpp
+++ b/cardmech/src/elasticity/assemble/DGElasticity.cpp
@@ -729,16 +729,30 @@ Vector DGElasticity::Projection(const Vector &u, const Vector &reference) const
}
auto r = projection.find_row(c());
int childCount = recursiveCells[dlevel].size();
- DGVectorFieldElement E(u, c);
- int nV = E.shape_size() * E.Dim();
- for (const auto &child: recursiveCells[dlevel]) {
- auto fineRow = reference.find_row(child->Center());
- for(int i=0;i<nV;++i) {
- projection(r, i) += reference(fineRow, i) / childCount;
+ DGVectorFieldElement E(projection, c);
+
+ for(int i = 0;i< E.NodalPoints(); ++i){
+ Point z = E.NodalPoint(i);
+ bool found=false;
+ for (const auto &child: recursiveCells[dlevel]) {
+ DGVectorFieldElement CE(reference, reference.find_cell(child->Center()));
+ for(int j = 0;j< CE.NodalPoints(); ++j){
+ if(CE.NodalPoint(j) == z){
+ auto rf = reference.find_row(child->Center());
+ for(int k=0;k<E.Dim();++k) {
+ projection(r, k*E.shape_size() + i) = reference(rf, k*CE.shape_size() + j);
+ }
+ found = true;
+ }
+ if(found)
+ break;
+ }
+ if(found)
+ break;
}
}
}
- return projection;
+ return projection;
}
@@ -815,7 +829,7 @@ double DGElasticity::EnergyError(const Vector &u) const {
double w = E.QWeight(q);
const Point &z = E.QPoint(q);
Tensor DD = (E.VectorGradient(q, u) - eProblem.DeformationGradient(timeStep, z));
- err += w * abs(cellMat.IsoDerivative(sym(One + DD), DD));
+ err += w * cellMat.IsoDerivative(sym(One + DD), DD);
err += w * abs(cellMat.AnisoDerivative(One + DD, Q, DD));
}
}
@@ -833,7 +847,7 @@ double DGElasticity::EnergyNorm(const Vector &u) const {
for (int q = 0; q < E.nQ(); ++q) {
double w = E.QWeight(q);
Tensor DD = E.VectorGradient(q, u);
- err += w * abs( cellMat.TotalDerivative(One + DD, One + DD, Q , DD));
+ err += w * abs(cellMat.TotalDerivative(One + DD, One+DD, Q, DD));
}
}
return sqrt(PPM->Sum(err));
diff --git a/cardmech/src/elasticity/problem/ElasticityBenchmarkProblem.hpp b/cardmech/src/elasticity/problem/ElasticityBenchmarkProblem.hpp
index 0ea525b56..47f205f54 100644
--- a/cardmech/src/elasticity/problem/ElasticityBenchmarkProblem.hpp
+++ b/cardmech/src/elasticity/problem/ElasticityBenchmarkProblem.hpp
@@ -50,7 +50,7 @@ public:
isPlanar(planar) {}
string MeshName() const override {
- return (isPlanar ? "UnitSquare" : "UnitCube");
+ return (isPlanar ? "UnitBlock2DQuad" : "UnitBlock3DQuad");
}
VectorField Deformation(const Point &z) const {
diff --git a/tools/jobsystem/run_dg_benchmarks.py b/tools/jobsystem/run_dg_benchmarks.py
new file mode 100644
index 000000000..e69de29bb
--
GitLab