From 5241c297127cf89ab64276eacfb8e6f4dcc1c405 Mon Sep 17 00:00:00 2001
From: acrovato <a.crovato@uliege.be>
Date: Wed, 16 Jun 2021 14:24:59 +0200
Subject: [PATCH] Remove Mem, mov methods to Element.
---
flow/src/wAdjoint.cpp | 2 +-
flow/src/wBlowingResidual.cpp | 4 +-
flow/src/wBoundary.cpp | 4 -
flow/src/wFreestreamResidual.cpp | 11 +-
flow/src/wKuttaElement.cpp | 3 -
flow/src/wKuttaResidual.cpp | 19 +--
flow/src/wLoadFunctional.cpp | 30 ++--
flow/src/wNewton.cpp | 5 +-
flow/src/wPotentialResidual.cpp | 43 +++---
flow/src/wProblem.cpp | 36 +++--
flow/src/wProblem.h | 4 +-
flow/src/wSolver.cpp | 2 +-
flow/src/wWakeElement.cpp | 3 -
flow/src/wWakeResidual.cpp | 39 +++---
heat/src/wHeatTerm.cpp | 18 +--
heat/src/wProblem.cpp | 13 +-
heat/src/wProblem.h | 2 +-
heat/src/wSolver.cpp | 4 +-
katoptron/src/LinearSolver.h | 2 +-
katoptron/src/ResultsDecl.hpp | 2 +-
katoptron/src/ResultsDef.hpp | 4 +-
katoptron/src/StressComputation.h | 12 +-
mirrors/src/wProblem.cpp | 15 +-
mirrors/src/wProblem.h | 2 +-
mirrors/src/wSolver.cpp | 6 +-
mirrors/src/wThermoMecaTerm.cpp | 22 ++-
tbox/_src/tboxw.h | 1 -
tbox/src/tbox.h | 6 -
tbox/src/wElement.cpp | 66 ++++++---
tbox/src/wElement.h | 66 +++++++--
tbox/src/wElement.inl | 166 ++++++++++++++++++++++
tbox/src/wHex8.cpp | 103 ++++++++++----
tbox/src/wHex8.h | 16 +--
tbox/src/wHex8.inl | 8 --
tbox/src/wLine2.cpp | 147 +++++++++++++-------
tbox/src/wLine2.h | 24 ++--
tbox/src/wLine2.inl | 8 --
tbox/src/wMem.cpp | 59 --------
tbox/src/wMem.h | 86 ------------
tbox/src/wMem.inl | 119 ----------------
tbox/src/wMemHex8.cpp | 68 ---------
tbox/src/wMemHex8.h | 50 -------
tbox/src/wMemLine2.cpp | 65 ---------
tbox/src/wMemLine2.h | 50 -------
tbox/src/wMemQuad4.cpp | 81 -----------
tbox/src/wMemQuad4.h | 50 -------
tbox/src/wMemTetra4.cpp | 68 ---------
tbox/src/wMemTetra4.h | 50 -------
tbox/src/wMemTri3.cpp | 81 -----------
tbox/src/wMemTri3.h | 50 -------
tbox/src/wMshCrack.cpp | 9 +-
tbox/src/wMshDeform.cpp | 22 ++-
tbox/src/wQuad4.cpp | 126 ++++++++++++-----
tbox/src/wQuad4.h | 20 ++-
tbox/src/wQuad4.inl | 8 --
tbox/src/wTetra4.cpp | 152 +++++++++++++-------
tbox/src/wTetra4.h | 20 ++-
tbox/src/wTetra4.inl | 8 --
tbox/src/wTri3.cpp | 223 +++++++++++++++++++-----------
tbox/src/wTri3.h | 28 ++--
tbox/src/wTri3.inl | 8 --
tbox/src/wUnitTest.cpp | 3 +-
waves/src/wProblem.cpp | 11 +-
waves/src/wProblem.h | 2 +-
waves/src/wTimeIntegration.cpp | 4 +-
waves/src/wWaveTerm.cpp | 9 +-
66 files changed, 971 insertions(+), 1477 deletions(-)
create mode 100644 tbox/src/wElement.inl
delete mode 100644 tbox/src/wMem.cpp
delete mode 100644 tbox/src/wMem.h
delete mode 100644 tbox/src/wMem.inl
delete mode 100644 tbox/src/wMemHex8.cpp
delete mode 100644 tbox/src/wMemHex8.h
delete mode 100644 tbox/src/wMemLine2.cpp
delete mode 100644 tbox/src/wMemLine2.h
delete mode 100644 tbox/src/wMemQuad4.cpp
delete mode 100644 tbox/src/wMemQuad4.h
delete mode 100644 tbox/src/wMemTetra4.cpp
delete mode 100644 tbox/src/wMemTetra4.h
delete mode 100644 tbox/src/wMemTri3.cpp
delete mode 100644 tbox/src/wMemTri3.h
diff --git a/flow/src/wAdjoint.cpp b/flow/src/wAdjoint.cpp
index 6774d1bb..a5d9be24 100644
--- a/flow/src/wAdjoint.cpp
+++ b/flow/src/wAdjoint.cpp
@@ -62,7 +62,7 @@ Adjoint::Adjoint(std::shared_ptr<Newton> _sol, std::shared_ptr<tbox::MshDeform>
verbose = sol->verbose;
// Update element gradients
- sol->pbl->initGradMem();
+ sol->pbl->initGradElems();
// Setup adjoint variables and gradients
lamClFlo.resize(sol->pbl->msh->nodes.size(), 0.);
diff --git a/flow/src/wBlowingResidual.cpp b/flow/src/wBlowingResidual.cpp
index 09dfcf30..15092468 100644
--- a/flow/src/wBlowingResidual.cpp
+++ b/flow/src/wBlowingResidual.cpp
@@ -20,7 +20,6 @@
#include "wElement.h"
#include "wCache.h"
#include "wGauss.h"
-#include "wMem.h"
using namespace tbox;
using namespace flow;
@@ -35,7 +34,6 @@ Eigen::VectorXd BlowingResidual::build(Element const &e, std::vector<double> con
// Get pre-computed values
Cache &cache = e.getVCache();
Gauss &gauss = cache.getVGauss();
- Mem &mem = e.getVMem();
// Blowing velocity
double vn = f.eval(e, phi, 0);
@@ -43,6 +41,6 @@ Eigen::VectorXd BlowingResidual::build(Element const &e, std::vector<double> con
// Integral of (normal) transpiration velocity
Eigen::VectorXd b = Eigen::VectorXd::Zero(e.nodes.size());
for (size_t k = 0; k < gauss.getN(); ++k)
- b += cache.getSf(k) * vn * gauss.getW(k) * mem.getDetJ(k);
+ b += cache.getSf(k) * vn * gauss.getW(k) * e.getDetJ(k);
return b;
}
diff --git a/flow/src/wBoundary.cpp b/flow/src/wBoundary.cpp
index 61e3bfb1..c01595b6 100644
--- a/flow/src/wBoundary.cpp
+++ b/flow/src/wBoundary.cpp
@@ -21,10 +21,6 @@
#include "wTag.h"
#include "wElement.h"
#include "wNode.h"
-#include "wFct0.h"
-#include "wCache.h"
-#include "wGauss.h"
-#include "wMem.h"
#include <unordered_set>
#include <iomanip>
#include <fstream>
diff --git a/flow/src/wFreestreamResidual.cpp b/flow/src/wFreestreamResidual.cpp
index afdd43ae..b2678b45 100644
--- a/flow/src/wFreestreamResidual.cpp
+++ b/flow/src/wFreestreamResidual.cpp
@@ -20,7 +20,6 @@
#include "wElement.h"
#include "wCache.h"
#include "wGauss.h"
-#include "wMem.h"
using namespace tbox;
using namespace flow;
@@ -35,15 +34,14 @@ Eigen::VectorXd FreestreamResidual::build(Element const &e, std::vector<double>
// Get pre-computed values
Cache &cache = e.getVCache();
Gauss &gauss = cache.getVGauss();
- Mem &mem = e.getVMem();
// Normal velocity
- double vn = fs.f.eval(e, phi, 0).dot(e.normal());
+ double vn = fs.f.eval(e, phi, 0).dot(e.getNormal());
// Integral of mass flux (density times normal velocity) through freestream boundary
Eigen::VectorXd b = Eigen::VectorXd::Zero(e.nodes.size());
for (size_t k = 0; k < gauss.getN(); ++k)
- b += cache.getSf(k) * vn * gauss.getW(k) * mem.getDetJ(k);
+ b += cache.getSf(k) * vn * gauss.getW(k) * e.getDetJ(k);
return b;
}
@@ -57,14 +55,13 @@ Eigen::VectorXd FreestreamResidual::buildGradientAoa(tbox::Element const &e, std
// Get pre-computed values
Cache &cache = e.getVCache();
Gauss &gauss = cache.getVGauss();
- Mem &mem = e.getVMem();
// Gradient of normal velocity
- double dvn = fs.f.evalD(e, phi, 0).dot(e.normal());
+ double dvn = fs.f.evalD(e, phi, 0).dot(e.getNormal());
// Build mass flux gradient
Eigen::VectorXd b = Eigen::VectorXd::Zero(e.nodes.size());
for (size_t k = 0; k < gauss.getN(); ++k)
- b += cache.getSf(k) * dvn * gauss.getW(k) * mem.getDetJ(k);
+ b += cache.getSf(k) * dvn * gauss.getW(k) * e.getDetJ(k);
return b;
}
diff --git a/flow/src/wKuttaElement.cpp b/flow/src/wKuttaElement.cpp
index 798fd778..c155c59f 100644
--- a/flow/src/wKuttaElement.cpp
+++ b/flow/src/wKuttaElement.cpp
@@ -16,9 +16,6 @@
#include "wKuttaElement.h"
#include "wNode.h"
-#include "wMem.h"
-#include "wCache.h"
-#include "wGauss.h"
using namespace tbox;
using namespace flow;
diff --git a/flow/src/wKuttaResidual.cpp b/flow/src/wKuttaResidual.cpp
index d7adcac9..ab56bd5e 100644
--- a/flow/src/wKuttaResidual.cpp
+++ b/flow/src/wKuttaResidual.cpp
@@ -20,7 +20,6 @@
#include "wElement.h"
#include "wCache.h"
#include "wGauss.h"
-#include "wMem.h"
using namespace tbox;
using namespace flow;
@@ -35,12 +34,11 @@ Eigen::MatrixXd KuttaResidual::buildFixed(KuttaElement const &ke, std::vector<do
// Get shape functions, Jacobian and Gauss points
Cache &sCache = ke.surE->getVCache();
Gauss &sGauss = sCache.getVGauss();
- Mem &sMem = ke.surE->getVMem();
// Get size
size_t nRow = ke.nRow;
// Velocity
- Eigen::RowVectorXd Aup = ke.volE->computeGradient(phi, 0).transpose() * ke.volE->getVMem().getJinv(0) * ke.volE->getVCache().getDsf(0);
+ Eigen::RowVectorXd Aup = ke.volE->computeGradient(phi, 0).transpose() * ke.volE->getJinv(0) * ke.volE->getVCache().getDsf(0);
// Build
Eigen::MatrixXd A = Eigen::MatrixXd::Zero(nRow, ke.nCol);
@@ -52,7 +50,7 @@ Eigen::MatrixXd KuttaResidual::buildFixed(KuttaElement const &ke, std::vector<do
for (size_t i = 0; i < nRow; ++i)
sfp(i) = sf(ke.teN[i].first);
// wake contribution
- A += ke.sign * sfp * Aup * sGauss.getW(k) * sMem.getDetJ(k);
+ A += ke.sign * sfp * Aup * sGauss.getW(k) * ke.surE->getDetJ(k);
}
return A;
}
@@ -67,7 +65,6 @@ Eigen::VectorXd KuttaResidual::build(KuttaElement const &ke, std::vector<double>
// Get shape functions, Jacobian and Gauss points
Cache &sCache = ke.surE->getVCache();
Gauss &sGauss = sCache.getVGauss();
- Mem &sMem = ke.surE->getVMem();
// Get size
size_t nRow = ke.nRow;
@@ -85,7 +82,7 @@ Eigen::VectorXd KuttaResidual::build(KuttaElement const &ke, std::vector<double>
for (size_t i = 0; i < nRow; ++i)
sfp(i) = sf(ke.teN[i].first);
// wake contribution
- b += ke.sign * sfp * dPhi2 * sGauss.getW(k) * sMem.getDetJ(k);
+ b += ke.sign * sfp * dPhi2 * sGauss.getW(k) * ke.surE->getDetJ(k);
}
return b;
}
@@ -115,10 +112,8 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> KuttaResidual::buildGradientMesh(Ku
Cache &surCache = ke.surE->getVCache();
Gauss &surGauss = surCache.getVGauss();
// Get Jacobian
- Mem &surMem = ke.surE->getVMem();
- Mem &volMem = ke.volE->getVMem();
- Eigen::MatrixXd const &iJ = volMem.getJinv(0);
- std::vector<Eigen::MatrixXd> const &dJ = volMem.getGradJ(0);
+ Eigen::MatrixXd const &iJ = ke.volE->getJinv(0);
+ std::vector<Eigen::MatrixXd> const &dJ = ke.volE->getGradJ(0);
// Get size
size_t nRow = ke.nRow;
size_t nDim = ke.nDim;
@@ -135,8 +130,8 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> KuttaResidual::buildGradientMesh(Ku
{
// weight and Jacobian determinant
double w = surGauss.getW(k);
- double detJ = surMem.getDetJ(k);
- std::vector<double> dDetJ = surMem.getGradDetJ(k);
+ double detJ = ke.surE->getDetJ(k);
+ std::vector<double> dDetJ = ke.surE->getGradDetJ(k);
// shape functions
Eigen::VectorXd const &sf = surCache.getSf(k);
Eigen::VectorXd sfp(nRow);
diff --git a/flow/src/wLoadFunctional.cpp b/flow/src/wLoadFunctional.cpp
index feefdaf9..6c0bbe17 100644
--- a/flow/src/wLoadFunctional.cpp
+++ b/flow/src/wLoadFunctional.cpp
@@ -20,7 +20,6 @@
#include "wElement.h"
#include "wCache.h"
#include "wGauss.h"
-#include "wMem.h"
using namespace tbox;
using namespace flow;
@@ -35,17 +34,16 @@ Eigen::VectorXd LoadFunctional::build(Element const &e, Element const &eV, std::
// Get pre-computed values
Cache &cache = e.getVCache();
Gauss &gauss = cache.getVGauss();
- Mem &mem = e.getVMem();
// Pressure force coefficient in associated volume
- Eigen::Vector3d fcp = cp.eval(eV, phi, 0) * e.normal();
+ Eigen::Vector3d fcp = cp.eval(eV, phi, 0) * e.getNormal();
// Integrate cp on the surface
Eigen::VectorXd b = Eigen::VectorXd::Zero(3 * e.nodes.size());
for (size_t k = 0; k < gauss.getN(); ++k)
{
Eigen::VectorXd sf = cache.getSf(k);
- double wdj = gauss.getW(k) * mem.getDetJ(k);
+ double wdj = gauss.getW(k) * e.getDetJ(k);
for (size_t i = 0; i < e.nodes.size(); ++i)
b.block<3, 1>(i * 3, 0) += sf(i) * fcp * wdj; // stack nodal force vectors [fx1, fy1, fz1, fx2, fy2, fz2, ...]
}
@@ -62,17 +60,16 @@ Eigen::MatrixXd LoadFunctional::buildGradientFlow(Element const &e, Element cons
// Get pre-computed values
Cache &cache = e.getVCache();
Gauss &gauss = cache.getVGauss();
- Mem &mem = e.getVMem();
// Gradient of pressure force coefficient in associated volume
- Eigen::MatrixXd dfcp = e.normal() * cp.evalD(eV, phi, 0) * eV.computeGradient(phi, 0).transpose() * eV.getVMem().getJinv(0) * eV.getVCache().getDsf(0);
+ Eigen::MatrixXd dfcp = e.getNormal() * cp.evalD(eV, phi, 0) * eV.computeGradient(phi, 0).transpose() * eV.getJinv(0) * eV.getVCache().getDsf(0);
// Integrate cp on the surface
Eigen::MatrixXd A = Eigen::MatrixXd::Zero(3 * e.nodes.size(), eV.nodes.size());
for (size_t k = 0; k < gauss.getN(); ++k)
{
Eigen::VectorXd sf = cache.getSf(k);
- double wdj = gauss.getW(k) * mem.getDetJ(k);
+ double wdj = gauss.getW(k) * e.getDetJ(k);
for (size_t i = 0; i < e.nodes.size(); ++i)
A.block(i * 3, 0, 3, eV.nodes.size()) += sf(i) * dfcp * wdj;
}
@@ -92,17 +89,16 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> LoadFunctional::buildGradientMesh(E
// Get pre-computed values
Cache &cache = e.getVCache();
Gauss &gauss = cache.getVGauss();
- Mem &mem = e.getVMem();
// Variables
- double cP = cp.eval(eV, phi, 0); // pressure coefficient
- Eigen::Vector3d nrm = e.normal(); // unit normal
+ double cP = cp.eval(eV, phi, 0); // pressure coefficient
+ Eigen::Vector3d nrm = e.getNormal(); // unit normal
// Gradients
- double dCp = cp.evalD(eV, phi, 0); // pressure coefficient
- Eigen::VectorXd dPhi = eV.computeGradient(phi, 0); // potential
- Eigen::MatrixXd const &iJ = eV.getVMem().getJinv(0); // inverse Jacobian
- std::vector<Eigen::MatrixXd> const &dJ = eV.getVMem().getGradJ(0); // Jacobian
- std::vector<Eigen::Vector3d> dNrm = e.gradientNormal(); // unit normal
+ double dCp = cp.evalD(eV, phi, 0); // pressure coefficient
+ Eigen::VectorXd dPhi = eV.computeGradient(phi, 0); // potential
+ Eigen::MatrixXd const &iJ = eV.getJinv(0); // inverse Jacobian
+ std::vector<Eigen::MatrixXd> const &dJ = eV.getGradJ(0); // Jacobian
+ std::vector<Eigen::Vector3d> dNrm = e.getGradNormal(); // unit normal
// Build element contributions
Eigen::MatrixXd As = Eigen::MatrixXd::Zero(3 * e.nodes.size(), e.nodes.size() * nDim);
@@ -111,8 +107,8 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> LoadFunctional::buildGradientMesh(E
{
Eigen::VectorXd sf = cache.getSf(k);
double w = gauss.getW(k);
- double detJ = mem.getDetJ(k);
- std::vector<double> dDetJ = mem.getGradDetJ(k);
+ double detJ = e.getDetJ(k);
+ std::vector<double> dDetJ = e.getGradDetJ(k);
for (size_t i = 0; i < e.nodes.size(); ++i)
{
for (size_t m = 0; m < nDim; ++m)
diff --git a/flow/src/wNewton.cpp b/flow/src/wNewton.cpp
index 71a50f75..ec0a579a 100644
--- a/flow/src/wNewton.cpp
+++ b/flow/src/wNewton.cpp
@@ -33,7 +33,6 @@
#include "wTag.h"
#include "wNode.h"
#include "wElement.h"
-#include "wMem.h"
#include "wLinesearch.h"
#include "wLinearSolver.h"
@@ -508,10 +507,10 @@ void Newton::findUpwd()
// -> find (cgU-cgE) \dot vel_ closest to one
double SPbest = -1;
size_t idx = 0;
- Eigen::VectorXd cgE = (p.first)->getVMem().getCg().block(0, 0, pbl->nDim, 1);
+ Eigen::VectorXd cgE = (p.first)->getCg().block(0, 0, pbl->nDim, 1);
for (size_t i = 0; i < p.second.size(); ++i)
{
- Eigen::VectorXd cgU = (p.second[i])->getVMem().getCg().block(0, 0, pbl->nDim, 1);
+ Eigen::VectorXd cgU = (p.second[i])->getCg().block(0, 0, pbl->nDim, 1);
double SP = (cgU - cgE).normalized().dot(vel_);
if (fabs(SP - 1) < fabs(SPbest - 1))
{
diff --git a/flow/src/wPotentialResidual.cpp b/flow/src/wPotentialResidual.cpp
index f98adbc5..1e683747 100644
--- a/flow/src/wPotentialResidual.cpp
+++ b/flow/src/wPotentialResidual.cpp
@@ -20,7 +20,6 @@
#include "wElement.h"
#include "wCache.h"
#include "wGauss.h"
-#include "wMem.h"
using namespace tbox;
using namespace flow;
@@ -35,16 +34,15 @@ Eigen::MatrixXd PotentialResidual::buildFixed(Element const &e, std::vector<doub
// Get pre-computed values
Cache &cache = e.getVCache();
Gauss &gauss = cache.getVGauss();
- Mem &mem = e.getVMem();
// Build
Eigen::MatrixXd A = Eigen::MatrixXd::Zero(e.nodes.size(), e.nodes.size());
for (size_t k = 0; k < gauss.getN(); ++k)
{
// Shape functions gradient
- Eigen::MatrixXd const &dSf = mem.getJinv(k) * cache.getDsf(k);
+ Eigen::MatrixXd const &dSf = e.getJinv(k) * cache.getDsf(k);
// Elementary stiffness matrix
- A += dSf.transpose() * dSf * fluid.rho.eval(e, phi, k) * gauss.getW(k) * mem.getDetJ(k);
+ A += dSf.transpose() * dSf * fluid.rho.eval(e, phi, k) * gauss.getW(k) * e.getDetJ(k);
}
return A;
}
@@ -60,12 +58,11 @@ Eigen::VectorXd PotentialResidual::build(Element const &e, Element const &eU, st
// Get pre-computed values
Cache &cache = e.getVCache();
Gauss &gauss = cache.getVGauss();
- Mem &mem = e.getVMem();
// Subsonic contribution
Eigen::VectorXd b = Eigen::VectorXd::Zero(e.nodes.size());
for (size_t k = 0; k < gauss.getN(); ++k)
- b += fluid.rho.eval(e, phi, k) * (mem.getJinv(k) * cache.getDsf(k)).transpose() * e.computeGradient(phi, k) * gauss.getW(k) * mem.getDetJ(k);
+ b += fluid.rho.eval(e, phi, k) * (e.getJinv(k) * cache.getDsf(k)).transpose() * e.computeGradient(phi, k) * gauss.getW(k) * e.getDetJ(k);
// Supersonic contribution
double mach = fluid.mach.eval(e, phi, 0);
@@ -77,7 +74,7 @@ Eigen::VectorXd PotentialResidual::build(Element const &e, Element const &eU, st
b *= 1 - mu;
// contribution of upwind element
for (size_t k = 0; k < gauss.getN(); ++k)
- b += mu * rhoU * (mem.getJinv(k) * cache.getDsf(k)).transpose() * e.computeGradient(phi, k) * gauss.getW(k) * mem.getDetJ(k);
+ b += mu * rhoU * (e.getJinv(k) * cache.getDsf(k)).transpose() * e.computeGradient(phi, k) * gauss.getW(k) * e.getDetJ(k);
}
return b;
}
@@ -97,16 +94,15 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientFlo
// Get pre-computed values
Cache &cache = e.getVCache();
Gauss &gauss = cache.getVGauss();
- Mem &mem = e.getVMem();
// Subsonic contribution
Eigen::MatrixXd A1 = Eigen::MatrixXd::Zero(e.nodes.size(), e.nodes.size());
for (size_t k = 0; k < gauss.getN(); ++k)
{
// Gauss point and determinant
- double wdj = gauss.getW(k) * mem.getDetJ(k);
+ double wdj = gauss.getW(k) * e.getDetJ(k);
// Shape functions and solution gradient
- Eigen::MatrixXd const &dSf = mem.getJinv(k) * cache.getDsf(k);
+ Eigen::MatrixXd const &dSf = e.getJinv(k) * cache.getDsf(k);
Eigen::VectorXd dPhi = e.computeGradient(phi, k);
// rho * grad_phi*grad_psi + drho * grad_phi*grad_psi
@@ -126,7 +122,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientFlo
Eigen::VectorXd dPhiU = eU.computeGradient(phi, 0);
double rhoU = fluid.rho.eval(eU, phi, 0);
double dRhoU = fluid.rho.evalD(eU, phi, 0);
- Eigen::MatrixXd const &dSfU = eU.getVMem().getJinv(0) * eU.getVCache().getDsf(0);
+ Eigen::MatrixXd const &dSfU = eU.getJinv(0) * eU.getVCache().getDsf(0);
// upwinding
A1 *= 1 - mu;
@@ -134,9 +130,9 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientFlo
for (size_t k = 0; k < gauss.getN(); ++k)
{
// Gauss point and determinant
- double wdj = gauss.getW(k) * mem.getDetJ(k);
+ double wdj = gauss.getW(k) * e.getDetJ(k);
// Shape functions and solution gradient
- Eigen::MatrixXd const &dSf = mem.getJinv(k) * cache.getDsf(k);
+ Eigen::MatrixXd const &dSf = e.getJinv(k) * cache.getDsf(k);
Eigen::VectorXd dPhi = e.computeGradient(phi, k);
// mu*drhoU*grad_phi*grad_psi
@@ -166,7 +162,6 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientMes
// Get pre-computed values
Cache &cache = e.getVCache();
Gauss &gauss = cache.getVGauss();
- Mem &mem = e.getVMem();
// Subsonic contribution
Eigen::MatrixXd A1 = Eigen::MatrixXd::Zero(e.nodes.size(), nDim * e.nodes.size());
@@ -177,13 +172,13 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientMes
double rho = fluid.rho.eval(e, phi, k);
double dRho = fluid.rho.evalD(e, phi, k);
Eigen::MatrixXd const &dSf = cache.getDsf(k);
- Eigen::MatrixXd const &iJ = mem.getJinv(k);
+ Eigen::MatrixXd const &iJ = e.getJinv(k);
// Jacobian gradients
- std::vector<Eigen::MatrixXd> const &dJ = mem.getGradJ(k);
- std::vector<double> const &dDetJ = mem.getGradDetJ(k);
+ std::vector<Eigen::MatrixXd> const &dJ = e.getGradJ(k);
+ std::vector<double> const &dDetJ = e.getGradDetJ(k);
// Gauss points and determinant
double w = gauss.getW(k);
- double detJ = mem.getDetJ(k);
+ double detJ = e.getDetJ(k);
for (size_t i = 0; i < e.nodes.size(); ++i)
{
for (size_t m = 0; m < nDim; ++m)
@@ -208,8 +203,8 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientMes
Eigen::VectorXd dPhiU = eU.computeGradient(phi, 0);
double rhoU = fluid.rho.eval(eU, phi, 0);
double dRhoU = fluid.rho.evalD(eU, phi, 0);
- Eigen::MatrixXd const &iJU = eU.getVMem().getJinv(0);
- std::vector<Eigen::MatrixXd> const &dJU = eU.getVMem().getGradJ(0);
+ Eigen::MatrixXd const &iJU = eU.getJinv(0);
+ std::vector<Eigen::MatrixXd> const &dJU = eU.getGradJ(0);
// upwinding
A1 *= 1 - mu;
@@ -221,13 +216,13 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientMes
double rho = fluid.rho.eval(e, phi, k);
double dM = fluid.mach.evalD(e, phi, k);
Eigen::MatrixXd const &dSf = cache.getDsf(k);
- Eigen::MatrixXd const &iJ = mem.getJinv(k);
+ Eigen::MatrixXd const &iJ = e.getJinv(k);
// Jacobian gradients
- std::vector<Eigen::MatrixXd> const &dJ = mem.getGradJ(k);
- std::vector<double> const &dDetJ = mem.getGradDetJ(k);
+ std::vector<Eigen::MatrixXd> const &dJ = e.getGradJ(k);
+ std::vector<double> const &dDetJ = e.getGradDetJ(k);
// Gauss points and determinant
double w = gauss.getW(k);
- double detJ = mem.getDetJ(k);
+ double detJ = e.getDetJ(k);
for (size_t m = 0; m < nDim; ++m)
{
for (size_t i = 0; i < e.nodes.size(); ++i)
diff --git a/flow/src/wProblem.cpp b/flow/src/wProblem.cpp
index 966064ca..7d3d069f 100644
--- a/flow/src/wProblem.cpp
+++ b/flow/src/wProblem.cpp
@@ -16,7 +16,6 @@
#include "wProblem.h"
#include "wElement.h"
-#include "wMem.h"
#include "wMedium.h"
#include "wBoundary.h"
#include "wAssign.h"
@@ -181,63 +180,62 @@ void Problem::add(std::shared_ptr<Blowing> b)
}
/**
- * @brief Initialize the elements memory (values)
- * @authors Adrien Crovato
+ * @brief Initialize the elements precomputed values
*/
-void Problem::initMem()
+void Problem::initElems()
{
// Update volume
for (auto e : medium->tag->elems)
- e->getVMem().update(true);
+ e->initValues(true);
// Update surface (boundaries)
for (auto surf : bnds)
for (auto e : surf->groups[0]->tag->elems)
- e->getVMem().update(false);
+ e->initValues(false);
// Update surface (Freestream B.C.)
for (auto surf : fBCs)
for (auto e : surf->tag->elems)
- e->getVMem().update(false);
+ e->initValues(false);
// Update surface (Wake B.C.)
for (auto surf : wBCs)
{
for (auto e : surf->groups[0]->tag->elems)
- e->getVMem().update(false);
+ e->initValues(false);
for (auto e : surf->groups[1]->tag->elems)
- e->getVMem().update(false);
+ e->initValues(false);
}
// Update surface (Blowing B.C.) if not already done
for (auto surf : bBCs)
for (auto e : surf->tag->elems)
- if (!(e->getVMem().isValid()))
- e->getVMem().update(false);
+ if (!(e->hasValues()))
+ e->initValues(false);
}
/**
- * @brief Initialize the elements memory (gradients)
+ * @brief Initialize the elements precomputed gradients
* @authors Adrien Crovato
*/
-void Problem::initGradMem()
+void Problem::initGradElems()
{
// Update volume
for (auto e : medium->tag->elems)
- e->getVMem().initGradients(true);
+ e->initGradients();
// Update surface (boundaries)
for (auto surf : bnds)
for (auto e : surf->groups[0]->tag->elems)
- e->getVMem().initGradients(false);
+ e->initGradients();
// Update surface (Wake B.C.)
for (auto surf : wBCs)
{
for (auto e : surf->groups[0]->tag->elems)
- e->getVMem().initGradients(false);
+ e->initGradients();
for (auto e : surf->groups[1]->tag->elems)
- e->getVMem().initGradients(false);
+ e->initGradients();
}
// Update surface (Blowing B.C.) if not already done
for (auto surf : bBCs)
for (auto e : surf->tag->elems)
- if (!(e->getVMem().isGradientValid()))
- e->getVMem().initGradients(false);
+ if (!(e->hasGradients()))
+ e->initGradients();
}
/**
diff --git a/flow/src/wProblem.h b/flow/src/wProblem.h
index 3f4b9b63..152aa1f2 100644
--- a/flow/src/wProblem.h
+++ b/flow/src/wProblem.h
@@ -78,8 +78,8 @@ public:
#ifndef SWIG
void check() const;
void pin();
- void initMem();
- void initGradMem();
+ void initElems();
+ void initGradElems();
virtual void write(std::ostream &out) const override;
#endif
};
diff --git a/flow/src/wSolver.cpp b/flow/src/wSolver.cpp
index d44e29dc..1b15b049 100644
--- a/flow/src/wSolver.cpp
+++ b/flow/src/wSolver.cpp
@@ -73,7 +73,7 @@ Solver::Solver(std::shared_ptr<Problem> _pbl, std::shared_ptr<LinearSolver> _lin
// Check the problem, pin a degree of freedom if necessary, and update element memory
pbl->check();
pbl->pin();
- pbl->initMem();
+ pbl->initElems();
// Assign a row (unknown index) for each node, and setup the local list
rows.resize(pbl->msh->nodes.size());
diff --git a/flow/src/wWakeElement.cpp b/flow/src/wWakeElement.cpp
index a75ccaa0..6146ed78 100644
--- a/flow/src/wWakeElement.cpp
+++ b/flow/src/wWakeElement.cpp
@@ -16,9 +16,6 @@
#include "wWakeElement.h"
#include "wNode.h"
-#include "wMem.h"
-#include "wCache.h"
-#include "wGauss.h"
using namespace tbox;
using namespace flow;
diff --git a/flow/src/wWakeResidual.cpp b/flow/src/wWakeResidual.cpp
index 0b111e23..12562c7c 100644
--- a/flow/src/wWakeResidual.cpp
+++ b/flow/src/wWakeResidual.cpp
@@ -20,7 +20,6 @@
#include "wElement.h"
#include "wCache.h"
#include "wGauss.h"
-#include "wMem.h"
using namespace tbox;
using namespace flow;
@@ -40,9 +39,8 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildFixed(WakeElemen
Eigen::MatrixXd const &vUpDsf = we.volUpE->getVCache().getDsf(0);
Eigen::MatrixXd const &vLwDsf = we.volLwE->getVCache().getDsf(0);
// Get Jacobian
- Mem &sMem = we.surUpE->getVMem();
- Eigen::MatrixXd const &vUpJ = we.volUpE->getVMem().getJinv(0);
- Eigen::MatrixXd const &vLwJ = we.volLwE->getVMem().getJinv(0);
+ Eigen::MatrixXd const &vUpJ = we.volUpE->getJinv(0);
+ Eigen::MatrixXd const &vLwJ = we.volLwE->getJinv(0);
// Get size
size_t nRow = we.nRow;
size_t nDim = we.nDim;
@@ -53,7 +51,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildFixed(WakeElemen
for (size_t j = 0; j < nDim; ++j)
dSf(i, j) = vUpDsf(j, we.vsMap.at(we.wkN[i].first));
Eigen::VectorXd dSfp(nRow);
- dSfp = 0.5 * sqrt(sMem.getVol()) * dSf * vUpJ.transpose() * vi;
+ dSfp = 0.5 * sqrt(we.surUpE->getVol()) * dSf * vUpJ.transpose() * vi;
// Velocity
Eigen::RowVectorXd Aup = we.volUpE->computeGradient(phi, 0).transpose() * vUpJ * vUpDsf;
Eigen::RowVectorXd Alw = we.volLwE->computeGradient(phi, 0).transpose() * vLwJ * vLwDsf;
@@ -69,8 +67,8 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildFixed(WakeElemen
for (size_t i = 0; i < nRow; ++i)
sfp(i) = dSfp(i) + sf(we.wkN[i].first);
// wake contribution
- Au += sfp * Aup * sGauss.getW(k) * sMem.getDetJ(k);
- Al += sfp * Alw * sGauss.getW(k) * sMem.getDetJ(k);
+ Au += sfp * Aup * sGauss.getW(k) * we.surUpE->getDetJ(k);
+ Al += sfp * Alw * sGauss.getW(k) * we.surUpE->getDetJ(k);
}
return std::make_tuple(Au, Al);
}
@@ -87,8 +85,6 @@ std::tuple<Eigen::VectorXd, Eigen::VectorXd> WakeResidual::build(WakeElement con
Cache &sCache = we.surUpE->getVCache();
Gauss &sGauss = sCache.getVGauss();
Eigen::MatrixXd const &vUpDsf = we.volUpE->getVCache().getDsf(0);
- // Get Jacobian
- Mem &sMem = we.surUpE->getVMem();
// Get size
size_t nRow = we.nRow;
size_t nDim = we.nDim;
@@ -99,7 +95,7 @@ std::tuple<Eigen::VectorXd, Eigen::VectorXd> WakeResidual::build(WakeElement con
for (size_t j = 0; j < nDim; ++j)
dSf(i, j) = vUpDsf(j, we.vsMap.at(we.wkN[i].first));
Eigen::VectorXd dSfp(nRow);
- dSfp = 0.5 * sqrt(sMem.getVol()) * dSf * we.volUpE->getVMem().getJinv(0).transpose() * vi;
+ dSfp = 0.5 * sqrt(we.surUpE->getVol()) * dSf * we.volUpE->getJinv(0).transpose() * vi;
// Velocity squared
Eigen::VectorXd dPhiU = we.volUpE->computeGradient(phi, 0);
Eigen::VectorXd dPhiL = we.volLwE->computeGradient(phi, 0);
@@ -117,8 +113,8 @@ std::tuple<Eigen::VectorXd, Eigen::VectorXd> WakeResidual::build(WakeElement con
for (size_t i = 0; i < nRow; ++i)
sfp(i) = dSfp(i) + sf(we.wkN[i].first);
// wake contribution
- bu += sfp * dPhi2u * sGauss.getW(k) * sMem.getDetJ(k);
- bl += sfp * dPhi2l * sGauss.getW(k) * sMem.getDetJ(k);
+ bu += sfp * dPhi2u * sGauss.getW(k) * we.surUpE->getDetJ(k);
+ bl += sfp * dPhi2l * sGauss.getW(k) * we.surUpE->getDetJ(k);
}
return std::make_tuple(bu, bl);
}
@@ -152,15 +148,12 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buil
Gauss &surGauss = sCache.getVGauss();
Eigen::MatrixXd const &dSfUp = we.volUpE->getVCache().getDsf(0);
// Get Jacobian
- Mem &sMem = we.surUpE->getVMem();
- Mem &vUpMem = we.volUpE->getVMem();
- Mem &vLwMem = we.volLwE->getVMem();
- Eigen::MatrixXd const &iJUp = vUpMem.getJinv(0);
- Eigen::MatrixXd const &iJLw = vLwMem.getJinv(0);
- std::vector<Eigen::MatrixXd> const &dJUp = vUpMem.getGradJ(0);
- std::vector<Eigen::MatrixXd> const &dJLw = vLwMem.getGradJ(0);
+ Eigen::MatrixXd const &iJUp = we.volUpE->getJinv(0);
+ Eigen::MatrixXd const &iJLw = we.volLwE->getJinv(0);
+ std::vector<Eigen::MatrixXd> const &dJUp = we.volUpE->getGradJ(0);
+ std::vector<Eigen::MatrixXd> const &dJLw = we.volLwE->getGradJ(0);
// Get surface gradient
- std::vector<double> const &dSurf = sMem.getGradVol();
+ std::vector<double> const &dSurf = we.surUpE->getGradVol();
// Get size
size_t nRow = we.nRow;
size_t nDim = we.nDim;
@@ -169,7 +162,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buil
size_t nSur = we.surUpE->nodes.size();
// Stabilization term
- double sqSurf = sqrt(sMem.getVol());
+ double sqSurf = sqrt(we.surUpE->getVol());
Eigen::MatrixXd dSf(nRow, nDim);
for (size_t i = 0; i < nRow; ++i)
for (size_t m = 0; m < nDim; ++m)
@@ -206,8 +199,8 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buil
{
// weight and Jacobian determinant
double w = surGauss.getW(k);
- double detJ = sMem.getDetJ(k);
- std::vector<double> dDetJ = sMem.getGradDetJ(k);
+ double detJ = we.surUpE->getDetJ(k);
+ std::vector<double> dDetJ = we.surUpE->getGradDetJ(k);
// shape functions and stabilization
Eigen::VectorXd const &sf = sCache.getSf(k);
Eigen::VectorXd sfp(nRow);
diff --git a/heat/src/wHeatTerm.cpp b/heat/src/wHeatTerm.cpp
index 057e139c..6d0238a6 100644
--- a/heat/src/wHeatTerm.cpp
+++ b/heat/src/wHeatTerm.cpp
@@ -19,7 +19,6 @@
#include "wElement.h"
#include "wCache.h"
#include "wGauss.h"
-#include "wMem.h"
#include "wFct0.h"
#include "wFct2.h"
@@ -47,17 +46,16 @@ Eigen::MatrixXd HeatTerm::build(Element const &e, std::vector<double> const &u,
else // true run - use MPI cached results
{
// Gauss integration
- Mem &mem = e.getVMem();
for (size_t k = 0; k < gauss.getN(); ++k)
{
// Jacobian inverse, shape functions and function evaluation
- Eigen::MatrixXd const &J = mem.getJinv(k);
+ Eigen::MatrixXd const &J = e.getJinv(k);
Eigen::MatrixXd const &dff = cache.getDsf(k);
Eigen::MatrixXd fk;
f.eval(e, u, k, fk, fake);
// Elementary stiffness matrix
- K += (fk * J * dff).transpose() * J * dff * gauss.getW(k) * mem.getDetJ(k);
+ K += (fk * J * dff).transpose() * J * dff * gauss.getW(k) * e.getDetJ(k);
}
}
return K;
@@ -71,12 +69,11 @@ Eigen::VectorXd HeatTerm::build(Element const &e, std::vector<double> const &u,
// Get precomputed values
Cache &cache = e.getVCache();
Gauss &gauss = cache.getVGauss();
- Mem &mem = e.getVMem();
// Gauss integration
Eigen::VectorXd s = Eigen::VectorXd::Zero(e.nodes.size());
for (size_t k = 0; k < gauss.getN(); ++k)
- s += cache.getSf(k) * f.eval(e, u, k) * gauss.getW(k) * mem.getDetJ(k);
+ s += cache.getSf(k) * f.eval(e, u, k) * gauss.getW(k) * e.getDetJ(k);
return s;
}
@@ -89,7 +86,6 @@ Eigen::VectorXd HeatTerm::build2(Element const &e, std::vector<double> const &u,
// Get precomputed values
Cache &cache = e.getVCache();
Gauss &gauss = cache.getVGauss();
- Mem &mem = e.getVMem();
// Gauss integration
Eigen::VectorXd qV = Eigen::VectorXd::Zero(e.nodes.size());
@@ -101,7 +97,7 @@ Eigen::VectorXd HeatTerm::build2(Element const &e, std::vector<double> const &u,
f.eval(e, u, k, fk, false);
// Elementary flux vector
- qV += (fk * e.computeGradient(u, k)).transpose() * mem.getJinv(k) * dff * gauss.getW(k) * mem.getDetJ(k);
+ qV += (fk * e.computeGradient(u, k)).transpose() * e.getJinv(k) * dff * gauss.getW(k) * e.getDetJ(k);
}
return qV;
}
@@ -114,7 +110,6 @@ Eigen::VectorXd HeatTerm::computeFlux(Element const &e, std::vector<double> cons
// Get precomputed values
Cache &cache = e.getVCache();
Gauss &gauss = cache.getVGauss();
- Mem &mem = e.getVMem();
// Gauss integration
Eigen::VectorXd qV = Eigen::Vector2d::Zero(); // 2 dimensions
@@ -125,7 +120,7 @@ Eigen::VectorXd HeatTerm::computeFlux(Element const &e, std::vector<double> cons
Eigen::MatrixXd fk(2, 2);
f.eval(e, u, k, fk, false);
- qV -= fk * gradk * gauss.getW(k) * mem.getDetJ(k);
+ qV -= fk * gradk * gauss.getW(k) * e.getDetJ(k);
}
return qV;
}
@@ -138,7 +133,6 @@ Eigen::MatrixXd HeatTerm::computeMatrix(Element const &e, std::vector<double> co
// Get precomputed values
Cache &cache = e.getVCache();
Gauss &gauss = cache.getVGauss();
- Mem &mem = e.getVMem();
// Gauss integration
Eigen::MatrixXd out = Eigen::Matrix2d::Zero(); // 2 dimensions, should match fk size
@@ -148,7 +142,7 @@ Eigen::MatrixXd HeatTerm::computeMatrix(Element const &e, std::vector<double> co
Eigen::MatrixXd fk;
f.eval(e, u, k, fk, false);
- out += fk * gauss.getW(k) * mem.getDetJ(k);
+ out += fk * gauss.getW(k) * e.getDetJ(k);
}
return out;
}
diff --git a/heat/src/wProblem.cpp b/heat/src/wProblem.cpp
index 997214de..717fd161 100644
--- a/heat/src/wProblem.cpp
+++ b/heat/src/wProblem.cpp
@@ -20,7 +20,7 @@
#include "wSource.h"
#include "wPeriodic.h" // sinon: warning C4150: deletion of pointer to incomplete type 'heat::Periodic'; no destructor called
#include "wElement.h"
-#include "wMem.h"
+
#include "wTag.h"
using namespace heat;
@@ -39,21 +39,20 @@ Problem::~Problem()
}
/**
- * @brief Update the elements memory (Jacobian)
- * @authors Adrien Crovato
+ * @brief Initialize the elements precomputed values
*/
-void Problem::updateMem()
+void Problem::initElems()
{
// Update volume Jacobian
for (auto vol : media)
for (auto e : vol->tag->elems)
- e->getVMem().update(true);
+ e->initValues(true);
// Update volume Jacobian (volume source)
for (auto vol : srcs)
for (auto e : vol->tag->elems)
- e->getVMem().update(true);
+ e->initValues(true);
// Update surface Jacobian (Neumann B.C.)
for (auto surf : bnds)
for (auto e : surf->tag->elems)
- e->getVMem().update(false);
+ e->initValues(false);
}
\ No newline at end of file
diff --git a/heat/src/wProblem.h b/heat/src/wProblem.h
index dd385f93..c961be28 100644
--- a/heat/src/wProblem.h
+++ b/heat/src/wProblem.h
@@ -57,7 +57,7 @@ public:
void add(std::shared_ptr<Periodic> p) { pdic = p; }
#ifndef SWIG
- void updateMem();
+ void initElems();
virtual void write(std::ostream &out) const override;
#endif
};
diff --git a/heat/src/wSolver.cpp b/heat/src/wSolver.cpp
index 685cabd1..35bab22c 100644
--- a/heat/src/wSolver.cpp
+++ b/heat/src/wSolver.cpp
@@ -62,7 +62,7 @@ Solver::Solver(std::shared_ptr<Problem> _pbl, std::shared_ptr<tbox::LinearSolver
save = true;
// update element memory
- pbl->updateMem();
+ pbl->initElems();
}
Solver::~Solver()
@@ -266,7 +266,7 @@ void Solver::start(MshExport *mshWriter)
// (flux moyen . volume) sur l'element
Eigen::VectorXd qV = HeatTerm::computeFlux(*e, T1, mat->k);
- double V = e->computeV();
+ double V = e->getVol();
if (save)
{
diff --git a/katoptron/src/LinearSolver.h b/katoptron/src/LinearSolver.h
index 7cc27fc8..0fe02a2f 100644
--- a/katoptron/src/LinearSolver.h
+++ b/katoptron/src/LinearSolver.h
@@ -8,7 +8,7 @@
#include "wTag.h"
#include "wMedium.h"
#include "wNode.h"
-#include "wMem.h"
+
#include "Teuchos_ParameterList.hpp"
#include "wProblem.h"
diff --git a/katoptron/src/ResultsDecl.hpp b/katoptron/src/ResultsDecl.hpp
index 4eeef4aa..9fac3532 100644
--- a/katoptron/src/ResultsDecl.hpp
+++ b/katoptron/src/ResultsDecl.hpp
@@ -8,7 +8,7 @@
#include "wTag.h"
#include "wMedium.h"
#include "wNode.h"
-#include "wMem.h"
+
#include "Teuchos_ParameterList.hpp"
#include "wProblem.h"
diff --git a/katoptron/src/ResultsDef.hpp b/katoptron/src/ResultsDef.hpp
index b4e5ab91..985cfa13 100644
--- a/katoptron/src/ResultsDef.hpp
+++ b/katoptron/src/ResultsDef.hpp
@@ -25,9 +25,9 @@ void writeResultsVTK(
for (auto e : pbl.msh->elems)
{
if (e->type() == tbox::ELTYPE::HEX8 || e->type() == tbox::ELTYPE::TETRA4) // 3D
- e->getVMem().update(true);
+ e->initValues(true);
//else // 2D
- // e->getVMem().update(false);
+ // e->initValues(false);
}
/* end temporary fix*/
diff --git a/katoptron/src/StressComputation.h b/katoptron/src/StressComputation.h
index 2870bd25..ec6c1e81 100644
--- a/katoptron/src/StressComputation.h
+++ b/katoptron/src/StressComputation.h
@@ -5,13 +5,11 @@
#include "wSfHex8.h"
#include "wGaussHex8.h"
#include "wCacheHex8.h"
-#include "wMemHex8.h"
#include "wTetra4.h"
#include "wSfTetra4.h"
#include "wGaussTetra4.h"
#include "wCacheTetra4.h"
-#include "wMemTetra4.h"
/**
* @brief Those functions compute the stress and strain at nodes and in average per element.
@@ -19,9 +17,9 @@
* such that it will be possible to use parallel for in the future.
*/
-void strain_at_gp(size_t n_n, tbox::Mem &mem, tbox::Cache &cache, size_t a, const Eigen::MatrixXd &u_e, Eigen::MatrixXd &Strain)
+void strain_at_gp(size_t n_n, tbox::Element &el, tbox::Cache &cache, size_t a, const Eigen::MatrixXd &u_e, Eigen::MatrixXd &Strain)
{
- Eigen::MatrixXd const &J = mem.getJinv(a);
+ Eigen::MatrixXd const &J = el.getJinv(a);
Eigen::MatrixXd const &dffi = cache.getDsf(a);
Eigen::Matrix3d du_d_Phi = Eigen::Matrix3d::Zero();
for (size_t i = 0; i < n_n; ++i)
@@ -57,7 +55,6 @@ void compute_stress(size_t n_n, const Eigen::MatrixXd &Strain, Eigen::MatrixXd &
void strain_stress_x_y_z(tbox::Element *elem, Eigen::MatrixXd &Strain, Eigen::MatrixXd &Stress, const Eigen::MatrixXd &H, const Eigen::MatrixXd &D, const std::vector<double> &x, const std::vector<double> &y, const std::vector<double> &z, const std::vector<double> &T, bool thermal = false)
{
tbox::Cache &cache = elem->getVCache();
- tbox::Mem &mem = elem->getVMem();
size_t n_gp = cache.getVGauss().getN();
size_t n_n = elem->nodes.size();
@@ -74,7 +71,7 @@ void strain_stress_x_y_z(tbox::Element *elem, Eigen::MatrixXd &Strain, Eigen::Ma
Strain = Eigen::Matrix3d::Zero();
for (size_t a = 0; a < n_gp; ++a)
{
- strain_at_gp(n_n, mem, cache, a, u_e, Strain);
+ strain_at_gp(n_n, *elem, cache, a, u_e, Strain);
}
Strain /= n_n;
@@ -84,7 +81,6 @@ void strain_stress_x_y_z(tbox::Element *elem, Eigen::MatrixXd &Strain, Eigen::Ma
void strain_stress_x_y_z_at_node(tbox::Element *elem, std::vector<Eigen::MatrixXd> &Strain_at_node, std::vector<Eigen::MatrixXd> &Stress_at_node, const Eigen::MatrixXd &H, const Eigen::MatrixXd &D, const std::vector<double> &x, const std::vector<double> &y, const std::vector<double> &z, const std::vector<double> &T, bool thermal = false)
{
tbox::Cache &cache = elem->getVCache();
- tbox::Mem &mem = elem->getVMem();
size_t n_gp = cache.getVGauss().getN();
size_t n_n = elem->nodes.size();
@@ -109,7 +105,7 @@ void strain_stress_x_y_z_at_node(tbox::Element *elem, std::vector<Eigen::MatrixX
for (size_t a = 0; a < n_gp; ++a)
{
Strain = Eigen::Matrix3d::Zero();
- strain_at_gp(n_n, mem, cache, a, u_e, Strain);
+ strain_at_gp(n_n, *elem, cache, a, u_e, Strain);
Strain_at_gp[a] = Strain;
compute_stress(n_n, Strain, Stress_at_gp[a], H, D, T, thermal);
}
diff --git a/mirrors/src/wProblem.cpp b/mirrors/src/wProblem.cpp
index a708feed..16bd8f13 100644
--- a/mirrors/src/wProblem.cpp
+++ b/mirrors/src/wProblem.cpp
@@ -23,7 +23,7 @@
#include "wTNeumann.h"
#include "wTSource.h"
#include "wElement.h"
-#include "wMem.h"
+
#include "wTag.h"
using namespace mirrors;
@@ -32,26 +32,25 @@ Problem::Problem(std::shared_ptr<MshData> _msh, std::string _problem_type, doubl
}
/**
- * @brief Update the elements memory (Jacobian)
- * @authors Adrien Crovato
+ * @brief Initialize the elements precomputed values
*/
-void Problem::updateMem()
+void Problem::initElems()
{
// Update volume Jacobian
for (auto vol : media)
for (auto e : vol->tag->elems)
- e->getVMem().update(true);
+ e->initValues(true);
// Update volume Jacobian (volumic sources)
for (auto vol : Tsrcs)
for (auto e : vol->tag->elems)
- e->getVMem().update(true);
+ e->initValues(true);
// Update surface Jacobian (Neumann B.C.)
for (auto surf : nuBCs)
for (auto e : surf->tag->elems)
- e->getVMem().update(false);
+ e->initValues(false);
for (auto surf : nTBCs)
for (auto e : surf->tag->elems)
- e->getVMem().update(false);
+ e->initValues(false);
}
void Problem::write(std::ostream &out) const
diff --git a/mirrors/src/wProblem.h b/mirrors/src/wProblem.h
index 21a61be6..0ce6de2b 100644
--- a/mirrors/src/wProblem.h
+++ b/mirrors/src/wProblem.h
@@ -62,7 +62,7 @@ public:
Problem(std::shared_ptr<MshData> _msh, std::string problem_type, double T_ref = 0);
#ifndef SWIG
- void updateMem();
+ void initElems();
virtual void write(std::ostream &out) const override;
#endif
};
diff --git a/mirrors/src/wSolver.cpp b/mirrors/src/wSolver.cpp
index ef21f944..7a01d5c4 100644
--- a/mirrors/src/wSolver.cpp
+++ b/mirrors/src/wSolver.cpp
@@ -50,7 +50,7 @@
#include <sstream>
#include <iostream>
/*
-// [RB] not needed anymore
+// [RB] not needed anymore
#ifdef max
#undef max // pour std::numeric_limits<std::streamsize>::max()
#endif
@@ -62,7 +62,7 @@ using namespace mirrors;
Solver::Solver(Problem &_pbl, std::shared_ptr<tbox::LinearSolver> _linsol, int _nthreads) : pbl(_pbl), linsol(_linsol), nthreads(_nthreads)
{
// Update element memory
- pbl.updateMem();
+ pbl.initElems();
}
void Solver::write(std::ostream &out) const
@@ -620,7 +620,7 @@ void Solver::start(MshExport *mshWriter)
Tag ¤t_tag = *(ms->tag);
- Eigen::Vector3d normal = current_tag.elems[0]->normal();
+ Eigen::Vector3d normal = current_tag.elems[0]->getNormal();
double max_distance = 0.;
diff --git a/mirrors/src/wThermoMecaTerm.cpp b/mirrors/src/wThermoMecaTerm.cpp
index d8a005e5..712542b2 100644
--- a/mirrors/src/wThermoMecaTerm.cpp
+++ b/mirrors/src/wThermoMecaTerm.cpp
@@ -19,7 +19,6 @@
#include "wElement.h"
#include "wCache.h"
#include "wGauss.h"
-#include "wMem.h"
#include "wNode.h"
using namespace mirrors;
@@ -33,7 +32,6 @@ Eigen::MatrixXd ThermoMecaTerm::buildM(Element const &e)
// Get precomputed values
Cache &cache = e.getVCache();
Gauss &gauss = cache.getVGauss();
- Mem &mem = e.getVMem();
// Mass matrix
Eigen::MatrixXd M = Eigen::MatrixXd::Zero(e.nodes.size(), e.nodes.size());
@@ -41,7 +39,7 @@ Eigen::MatrixXd ThermoMecaTerm::buildM(Element const &e)
{
// Shape functions
Eigen::VectorXd const &ff = cache.getSf(k);
- M += ff * ff.transpose() * gauss.getW(k) * mem.getDetJ(k);
+ M += ff * ff.transpose() * gauss.getW(k) * e.getDetJ(k);
}
return M;
}
@@ -54,7 +52,6 @@ Eigen::MatrixXd ThermoMecaTerm::buildK(tbox::Element const &e, Eigen::MatrixXd c
// Get precomputed values
Cache &cache = e.getVCache();
Gauss &gauss = cache.getVGauss();
- Mem &mem = e.getVMem();
size_t ns = e.nodes.size(); // number of nodes
// Meca matrix
@@ -62,7 +59,7 @@ Eigen::MatrixXd ThermoMecaTerm::buildK(tbox::Element const &e, Eigen::MatrixXd c
for (size_t k = 0; k < gauss.getN(); ++k)
{
// Jacobian inverse
- Eigen::Matrix3d const &invJ = mem.getJinv(k);
+ Eigen::Matrix3d const &invJ = e.getJinv(k);
// Fill B matrix (shape functions)
Eigen::MatrixXd const &dff = cache.getDsf(k);
Eigen::MatrixXd B = Eigen::MatrixXd::Zero(6, 3 * ns);
@@ -82,7 +79,7 @@ Eigen::MatrixXd ThermoMecaTerm::buildK(tbox::Element const &e, Eigen::MatrixXd c
}
// Compute stiffness matrix
- K += B.transpose() * H * B * gauss.getW(k) * mem.getDetJ(k);
+ K += B.transpose() * H * B * gauss.getW(k) * e.getDetJ(k);
}
return K;
}
@@ -95,15 +92,14 @@ Eigen::MatrixXd ThermoMecaTerm::buildL(Element const &e, Eigen::MatrixXd const &
// Get precomputed values
Cache &cache = e.getVCache();
Gauss &gauss = cache.getVGauss();
- Mem &mem = e.getVMem();
// Thermal matrix
Eigen::MatrixXd L = Eigen::MatrixXd::Zero(e.nodes.size(), e.nodes.size());
for (size_t k = 0; k < gauss.getN(); ++k)
{
// Jacobian
- double detJ = mem.getDetJ(k);
- Eigen::Matrix3d const &J = mem.getJinv(k);
+ double detJ = e.getDetJ(k);
+ Eigen::Matrix3d const &J = e.getJinv(k);
Eigen::MatrixXd const &dff = cache.getDsf(k);
// [AC] I refactored this from an entry by entry assembly to a full matrix multiplication,
@@ -121,7 +117,6 @@ Eigen::MatrixXd ThermoMecaTerm::buildS(Element const &e, Eigen::MatrixXd const &
// Get precomputed values
Cache &cache = e.getVCache();
Gauss &gauss = cache.getVGauss();
- Mem &mem = e.getVMem();
size_t ns = e.nodes.size(); // number of nodes
// TM matrix
@@ -129,8 +124,8 @@ Eigen::MatrixXd ThermoMecaTerm::buildS(Element const &e, Eigen::MatrixXd const &
for (size_t k = 0; k < gauss.getN(); ++k)
{
// Jacobian and shape functions
- double detJ = mem.getDetJ(k);
- Eigen::Matrix3d const &J = mem.getJinv(k);
+ double detJ = e.getDetJ(k);
+ Eigen::Matrix3d const &J = e.getJinv(k);
Eigen::MatrixXd const &dff = cache.getDsf(k);
// S
@@ -148,7 +143,6 @@ void ThermoMecaTerm::strain_stress(Element const &e, Eigen::MatrixXd const &H, s
// Get precomputed values
Cache &cache = e.getVCache();
Gauss &gauss = cache.getVGauss();
- Mem &mem = e.getVMem();
size_t ns = e.nodes.size(); // number of nodes
Strain = Stress = Eigen::Matrix3d::Zero(); // 3 dimensions
@@ -163,7 +157,7 @@ void ThermoMecaTerm::strain_stress(Element const &e, Eigen::MatrixXd const &H, s
for (size_t k = 0; k < gauss.getN(); ++k)
{
// Jacobian and strain
- Eigen::Matrix3d const &J = mem.getJinv(k);
+ Eigen::Matrix3d const &J = e.getJinv(k);
Eigen::MatrixXd const &dff = cache.getDsf(k);
Eigen::Matrix3d du_d_Phi = Eigen::Matrix3d::Zero();
for (size_t i = 0; i < ns; ++i)
diff --git a/tbox/_src/tboxw.h b/tbox/_src/tboxw.h
index 85f85289..b9feacda 100644
--- a/tbox/_src/tboxw.h
+++ b/tbox/_src/tboxw.h
@@ -43,7 +43,6 @@
#include "wLazy.h"
#include "wLine2.h"
#include "wLinearSolver.h"
-#include "wMemTri3.h"
#include "wMshConvert.h"
#include "wMshCrack.h"
#include "wMshData.h"
diff --git a/tbox/src/tbox.h b/tbox/src/tbox.h
index 4c2304c6..b6814d4c 100644
--- a/tbox/src/tbox.h
+++ b/tbox/src/tbox.h
@@ -44,7 +44,6 @@ class Tag;
// Element handling
class Cache;
-class Mem;
// Point1
class Point1;
@@ -54,35 +53,30 @@ class CachePoint1;
class SfLine2;
class CacheLine2;
class GaussLine2;
-class MemLine2;
class Line2;
// Tri3
class SfTri3;
class CacheTri3;
class GaussTri3;
-class MemTri3;
class Tri3;
// Quad4
class SfQuad4;
class CacheQuad4;
class GaussQuad4;
-class MemQuad4;
class Quad4;
// Tetra4
class SfTetra4;
class CacheTetra4;
class GaussTetra4;
-class MemTetra4;
class Tetra4;
// Hex8
class SfHex8;
class CacheHex8;
class GaussHex8;
-class MemHex8;
class Hex8;
// Mesh utilities
diff --git a/tbox/src/wElement.cpp b/tbox/src/wElement.cpp
index aa8fba02..2aab2292 100644
--- a/tbox/src/wElement.cpp
+++ b/tbox/src/wElement.cpp
@@ -17,8 +17,6 @@
#include "wElement.h"
#include "wNode.h"
#include "wTag.h"
-#include "wFct0.h"
-#include "wMem.h"
#include "wCache.h"
#include "std_extra.h"
@@ -61,8 +59,8 @@ operator<<(std::ostream &out, ELTYPE const &obj)
Element::Element(int n, Tag *_ptag, Tag *_etag,
std::vector<Tag *> &_parts, std::vector<Node *> &nods)
- : wObject(), order(2), no(n), ptag(_ptag), etag(_etag), parts(_parts),
- nodes(nods)
+ : wObject(), order(0), hasVals(false), hasGrads(false), isVol(false),
+ no(n), ptag(_ptag), etag(_etag), parts(_parts), nodes(nods)
{
if (ptag)
ptag->elems.push_back(this);
@@ -73,37 +71,69 @@ Element::Element(int n, Tag *_ptag, Tag *_etag,
p->elems.push_back(this);
}
-Mem &Element::getVMem() const
+void Element::initValues(bool isvol, int ordr)
{
- throw std::runtime_error("Element::getVMem not implemented\n");
+ throw std::runtime_error("Element::initValues() not implemented!\n");
+}
+void Element::initGradients()
+{
+ throw std::runtime_error("Element::initGradients() not implemented!\n");
+}
+void Element::update()
+{
+ throw std::runtime_error("Element::update() not implemented!\n");
}
-double Element::buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const
+/**
+ * @brief Check that determinant of Jacobian is positive
+ */
+void Element::checkJac(size_t k) const
{
- throw std::runtime_error("Element::buildJ not implemented\n");
+ if (detJs[k] <= 0)
+ {
+ std::stringstream err;
+ err << "Negative jacobian detected for element: " << *this << "\n";
+ throw std::runtime_error(err.str());
+ }
}
-double Element::buildJ(size_t k) const
+
+void Element::buildJ()
{
throw std::runtime_error("Element::buildJ not implemented\n");
}
-void Element::buildGradientJ(size_t k, std::vector<double> &dDetJ, std::vector<Eigen::MatrixXd> &dJ) const
+void Element::buildSurfaceJ()
{
- throw std::runtime_error("Element::buildGradientJ not implemented\n");
+ throw std::runtime_error("Element::buildSurfaceJ not implemented\n");
}
-void Element::buildGradientJ(size_t k, std::vector<double> &dDetJ) const
+void Element::buildGradientJ()
{
throw std::runtime_error("Element::buildGradientJ not implemented\n");
}
+void Element::buildGradientSurfaceJ()
+{
+ throw std::runtime_error("Element::buildGradientSurfaceJ not implemented\n");
+}
-double Element::computeV() const
+void Element::computeV()
{
throw std::runtime_error("Element::computeV not implemented\n");
}
-std::vector<double> Element::buildGradientV(int dim) const
+void Element::buildGradientV(int dim)
{
throw std::runtime_error("Element::buildGradientV not implemented\n");
}
+/**
+ * @brief Compute center of gravity of element
+ */
+void Element::computeCg()
+{
+ cg = Eigen::Vector3d::Zero();
+ for (auto n : nodes)
+ cg += n->pos;
+ cg /= nodes.size();
+}
+
std::vector<Eigen::Vector3d> Element::computeTangents() const
{
throw std::runtime_error("Element::computeTangents not implemented\n");
@@ -113,13 +143,13 @@ std::vector<Eigen::Vector3d> Element::buildTangents(size_t k) const
throw std::runtime_error("Element::buildTangents not implemented\n");
}
-Eigen::Vector3d Element::normal() const
+void Element::computeNormal()
{
- throw std::runtime_error("Element::normal not implemented\n");
+ throw std::runtime_error("Element::computeNormal not implemented\n");
}
-std::vector<Eigen::Vector3d> Element::gradientNormal() const
+void Element::computeGradientNormal()
{
- throw std::runtime_error("Element::gradientNormal not implemented\n");
+ throw std::runtime_error("Element::computeGradientNormal not implemented\n");
}
/**
diff --git a/tbox/src/wElement.h b/tbox/src/wElement.h
index 72c4dd9b..3c02267a 100644
--- a/tbox/src/wElement.h
+++ b/tbox/src/wElement.h
@@ -49,8 +49,27 @@ TBOX_API std::ostream &operator<<(std::ostream &out, ELTYPE const &obj);
*/
class TBOX_API Element : public fwk::wObject
{
+ friend class MshCrack; // so that MshCrack can access normal and cg before update
protected:
int order; ///< order of integration (Gauss quadrature)
+ // Flags
+ bool hasVals; ///< true if (precomputed) values have been initialized
+ bool hasGrads; ///< true if (precomputed) gradients have been initialized
+ bool isVol; ///< true if element is a volume element
+ // Jacobian
+ std::vector<double> detJs; ///< determinant of Jacobian matrix
+ std::vector<Eigen::MatrixXd> Js; ///< Jacobian matrix
+ std::vector<Eigen::MatrixXd> iJs; ///< inverse of Jacobian matrix
+ std::vector<std::vector<double>> gradDetJs; ///< gradient of Jacobian determinant wrt to nodes position
+ std::vector<std::vector<Eigen::MatrixXd>> gradJs; ///< gradient of Jacobian matrix wrt to nodes position
+ // Volume and CG
+ double vol; ///< element volume
+ Eigen::Vector3d cg; ///< element center of gravity
+ std::vector<double> gradVol; ///< gradient of volume wrt to nodes position
+ // Normal vectors
+ Eigen::Vector3d normal; ///< unit normal vector
+ std::vector<Eigen::Vector3d> gradNormal; ///< graident of unit normal vector wrt to nodes position
+
public:
int no; ///< label
Tag *ptag; ///< physical tag (group)
@@ -70,30 +89,49 @@ public:
#ifndef SWIG
protected:
// Jacobian
- virtual double buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const;
- virtual double buildJ(size_t k) const;
- virtual void buildGradientJ(size_t k, std::vector<double> &dDetJ, std::vector<Eigen::MatrixXd> &dJ) const;
- virtual void buildGradientJ(size_t k, std::vector<double> &dDetJ) const;
+ void checkJac(size_t k) const;
+ virtual void buildJ();
+ virtual void buildSurfaceJ();
+ virtual void buildGradientJ();
+ virtual void buildGradientSurfaceJ();
+ // CG and Volume
+ void computeCg();
+ virtual void computeV();
+ virtual void buildGradientV(int dim);
+ // Tangent and normal vectors
+ virtual std::vector<Eigen::Vector3d> computeTangents() const;
+ virtual std::vector<Eigen::Vector3d> buildTangents(size_t k) const;
+ virtual void computeNormal();
+ virtual void computeGradientNormal();
public:
+ // Initialize and update precomputed values
+ virtual void initValues(bool isvol, int ordr = 2);
+ virtual void initGradients();
+ virtual void update();
// Getters
virtual Cache &getVCache() const = 0;
- virtual Mem &getVMem() const;
- // Volume (@todo private)
- virtual double computeV() const;
- virtual std::vector<double> buildGradientV(int dim) const;
- // Tangent and normal vectors (@todo private)
- virtual std::vector<Eigen::Vector3d> computeTangents() const;
- virtual std::vector<Eigen::Vector3d> buildTangents(size_t k) const;
- virtual Eigen::Vector3d normal() const;
- virtual std::vector<Eigen::Vector3d> gradientNormal() const;
+ inline bool hasValues() const;
+ inline bool hasGradients() const;
+ inline double getDetJ(size_t k) const;
+ inline Eigen::MatrixXd const &getJ(size_t k) const;
+ inline Eigen::MatrixXd const &getJinv(size_t k) const;
+ inline std::vector<double> const &getGradDetJ(size_t k) const;
+ inline std::vector<Eigen::MatrixXd> const &getGradJ(size_t k) const;
+ inline double getVol() const;
+ inline Eigen::Vector3d const &getCg() const;
+ inline std::vector<double> const &getGradVol() const;
+ inline Eigen::Vector3d const &getNormal() const;
+ inline std::vector<Eigen::Vector3d> const &getGradNormal() const;
// Utilities
virtual Eigen::VectorXd computeGradient(std::vector<double> const &u, size_t k) const;
- virtual Eigen::MatrixXd gp2node(Eigen::MatrixXd &values_at_gp) const;
+ Eigen::MatrixXd gp2node(Eigen::MatrixXd &values_at_gp) const;
virtual void write(std::ostream &out) const override;
#endif
};
+#include "wElement.inl"
+
} // namespace tbox
#endif //WELEMENT_H
diff --git a/tbox/src/wElement.inl b/tbox/src/wElement.inl
new file mode 100644
index 00000000..9bed1207
--- /dev/null
+++ b/tbox/src/wElement.inl
@@ -0,0 +1,166 @@
+/*
+ * Copyright 2020 University of Liège
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+/**
+ * @brief Return wether the element values have been initialized or not
+ */
+inline bool Element::hasValues() const
+{
+ return hasVals;
+}
+
+/**
+ * @brief Return wether the element gradients have been initialized or not
+ */
+inline bool Element::hasGradients() const
+{
+ return hasGrads;
+}
+
+/**
+ * @brief Return the determinant of the Jacobian at Gauss point k
+ */
+inline double Element::getDetJ(size_t k) const
+{
+ if (hasVals)
+ return detJs[k];
+ else
+ throw std::runtime_error("Element::getDetJ is invalid because element memory is not up-to-date!\n");
+}
+
+/**
+ * @brief Return the the Jacobian matrix at Gauss point k
+ */
+inline Eigen::MatrixXd const &Element::getJ(size_t k) const
+{
+ if (hasVals)
+ {
+ if (isVol)
+ return Js[k];
+ else
+ throw std::runtime_error("Element::getJ Jacobian matrix is not available for surface element!\n");
+ }
+ else
+ throw std::runtime_error("Element::getJ is invalid because element memory is not up-to-date!\n");
+}
+
+/**
+ * @brief Return the inverse of the Jacobian matrix at Gauss point k
+ */
+inline Eigen::MatrixXd const &Element::getJinv(size_t k) const
+{
+ if (hasVals)
+ {
+ if (isVol)
+ return iJs[k];
+ else
+ throw std::runtime_error("Element::getJinv inverse Jacobian matrix is not available for surface element!\n");
+ }
+ else
+ throw std::runtime_error("Element::getJinv is invalid because element memory is not up-to-date!\n");
+}
+
+/**
+ * @brief Return the volume of the element
+ */
+inline double Element::getVol() const
+{
+ if (hasVals)
+ return vol;
+ else
+ throw std::runtime_error("Mem::getVol is invalid because element memory is not up-to-date!\n");
+}
+
+/**
+ * @brief Return the center of gravity of the element
+ */
+inline Eigen::Vector3d const &Element::getCg() const
+{
+ if (hasVals)
+ return cg;
+ else
+ throw std::runtime_error("Element::getCg is invalid because element memory is not up-to-date!\n");
+}
+
+/**
+ * @brief Return the unit normal vector of the element
+ */
+inline Eigen::Vector3d const &Element::getNormal() const
+{
+ if (hasVals)
+ {
+ if (!isVol)
+ return normal;
+ else
+ throw std::runtime_error("Element::getNormal normal vector is not available for volume element!\n");
+ }
+ else
+ throw std::runtime_error("Element::getNormal is invalid because element memory is not up-to-date!\n");
+}
+
+/**
+ * @brief Return the gradient of determinant of the Jacobian with respect to node positions at Gauss point k
+ */
+inline std::vector<double> const &Element::getGradDetJ(size_t k) const
+{
+ if (hasGrads)
+ return gradDetJs[k];
+ else
+ throw std::runtime_error("Element::getGradDetJ is invalid because element memory (gradient) is not up-to-date!\n");
+}
+
+/**
+ * @brief Return the gradient of the Jacobian matrix with respect to node positions at Gauss point k
+ */
+inline std::vector<Eigen::MatrixXd> const &Element::getGradJ(size_t k) const
+{
+ if (hasGrads)
+ {
+ if (isVol)
+ return gradJs[k];
+ else
+ throw std::runtime_error("Element::getGradJ gradient of Jacobian matrix is not available for surface element!\n");
+ }
+ else
+ throw std::runtime_error("Element::getGradJ is invalid because element memory (gradient) is not up-to-date!\n");
+}
+
+/**
+ * @brief Return the gradient of the volume with respect to node positions
+ */
+inline std::vector<double> const &Element::getGradVol() const
+{
+ if (hasGrads)
+ return gradVol;
+ else
+ throw std::runtime_error("Element::getGradVol is invalid because element memory (gradient) is not up-to-date!\n");
+}
+
+/**
+ * @brief Return the gradient of the unit normal vector with respect to node positions
+ */
+inline std::vector<Eigen::Vector3d> const &Element::getGradNormal() const
+{
+ if (hasGrads)
+ {
+ if (!isVol)
+ return gradNormal;
+ else
+ throw std::runtime_error("Element::getGradNormal normal vector gradient is not available for volume element!\n");
+ }
+ else
+ throw std::runtime_error("Element::getGradNormal is invalid because element memory (gradient) is not up-to-date!\n");
+}
diff --git a/tbox/src/wHex8.cpp b/tbox/src/wHex8.cpp
index ae4c523a..606ce106 100644
--- a/tbox/src/wHex8.cpp
+++ b/tbox/src/wHex8.cpp
@@ -18,7 +18,6 @@
#include "wSfHex8.h"
#include "wGaussHex8.h"
#include "wCacheHex8.h"
-#include "wMemHex8.h"
#include "wNode.h"
using namespace tbox;
@@ -26,60 +25,108 @@ using namespace tbox;
Hex8::Hex8(int n, Tag *_ptag, Tag *_etag,
std::vector<Tag *> &_parts, std::vector<Node *> &nods) : Element(n, _ptag, _etag, _parts, nods)
{
- pmem = new MemHex8(*this, static_cast<int>(getCache(order).gauss.getN()));
}
-Hex8::~Hex8()
+
+/**
+ * @brief Return cache handling shape functions and Gauss points for all Hex8 elements
+ */
+Cache &Hex8::getVCache() const
{
- if (pmem)
- delete pmem;
+ return getCache(order);
}
/**
- * @brief Return private memory handling Jacobian for this Hex8 element
+ * @brief Initialize precomputed values
*/
-Mem &Hex8::getVMem() const
+void Hex8::initValues(bool isvol, int ordr)
{
- return getMem();
+ // Set flags
+ isVol = isvol;
+ order = ordr;
+ hasVals = true;
+ // Update values
+ size_t ngp = getCache(order).getVGauss().getN();
+ Js.resize(ngp);
+ iJs.resize(ngp);
+ detJs.resize(ngp);
+ this->update();
}
-
/**
- * @brief Return cache handling shape functions and Gauss points for all Hex8 elements
+ * @brief Initialize precomputed gradients
*/
-Cache &Hex8::getVCache() const
+void Hex8::initGradients()
{
- return getCache(order);
+ // Set flags
+ hasGrads = true;
+ // Update values
+ size_t ngp = getCache(order).getVGauss().getN();
+ gradJs.resize(ngp);
+ gradDetJs.resize(ngp);
+ this->update();
+}
+/**
+ * @brief Update Hex8 precomputed values and gradients
+ */
+void Hex8::update()
+{
+ // Update values
+ if (hasVals)
+ {
+ // Jacobian
+ if (!isVol)
+ throw std::runtime_error("Hex8::update(): the problem must always have the same dimension as Hex8!\n");
+ // Hex8 is always of same dimension as Problem (volume element)
+ this->buildJ();
+ // Volume and CG
+ this->computeV();
+ this->computeCg();
+ }
+ // Update gradients
+ if (hasGrads)
+ {
+ // Jacobian
+ if (!isVol)
+ throw std::runtime_error("Hex8::updateGradients(): the problem must always have the same dimension as Hex8!\n");
+ // Hex8 is always of same dimension as Problem (volume element)
+ this->buildGradientJ();
+ // Volume
+ this->buildGradientV(3);
+ }
}
/**
- * @brief Compute Jacobian matrix, inverse and determinant at Gauss point k
- *
+ * @brief Compute Jacobian matrix, inverse and determinant
+ *
* J_ij(3,3) = diN_n * xj_n (n is the node)
*/
-double Hex8::buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const
+void Hex8::buildJ()
{
- Eigen::MatrixXd const &dff = getCache(order).getDsf(k);
+ CacheHex8 &cache = getCache(order);
+ for (size_t k = 0; k < cache.getVGauss().getN(); ++k)
+ {
+ Eigen::MatrixXd const &dff = cache.getDsf(k);
- Eigen::Matrix3d JJ = Eigen::Matrix3d::Zero(); // temporary fixed-size matrix to efficiently compute the inverse
- size_t i = 0;
- for (auto it = nodes.begin(); it != nodes.end(); ++it, ++i)
- JJ += dff.col(i) * (*it)->pos.transpose();
- J = JJ;
- iJ = JJ.inverse();
- return JJ.determinant();
+ Eigen::Matrix3d JJ = Eigen::Matrix3d::Zero(); // temporary fixed-size matrix to efficiently compute the inverse
+ size_t i = 0;
+ for (auto it = nodes.begin(); it != nodes.end(); ++it, ++i)
+ JJ += dff.col(i) * (*it)->pos.transpose();
+ Js[k] = JJ;
+ iJs[k] = JJ.inverse();
+ detJs[k] = JJ.determinant();
+ this->checkJac(k);
+ }
}
/**
* @brief Compute the element volume
*/
-double Hex8::computeV() const
+void Hex8::computeV()
{
CacheHex8 &cache = getCache(order);
GaussHex8 &gauss = cache.gauss;
- MemHex8 &mem = getMem();
// Gauss integration
- double V = 0.0;
+ vol = 0.0;
for (size_t k = 0; k < gauss.getN(); ++k)
- V += gauss.getW(k) * mem.getDetJ(k);
- return V;
+ vol += gauss.getW(k) * this->getDetJ(k);
}
diff --git a/tbox/src/wHex8.h b/tbox/src/wHex8.h
index bf39aaa5..391b73dc 100644
--- a/tbox/src/wHex8.h
+++ b/tbox/src/wHex8.h
@@ -32,28 +32,24 @@ namespace tbox
class TBOX_API Hex8 : public Element
{
#ifndef SWIG
- friend class MemHex8;
-
private:
- MemHex8 *pmem; ///< private memory of precalculated values
inline static CacheHex8 &getCache(int n);
- inline MemHex8 &getMem() const;
protected:
- virtual double buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const override;
-
-public: // @todo
- virtual double computeV() const override;
+ virtual void buildJ() override;
+ virtual void computeV() override;
#endif
public:
Hex8(int n, Tag *_ptag, Tag *_etag, std::vector<Tag *> &_parts, std::vector<Node *> &nods);
- virtual ~Hex8();
+ virtual ~Hex8() {}
virtual ELTYPE type() const override { return ELTYPE::HEX8; }
#ifndef SWIG
+ virtual void initValues(bool isvol, int ordr = 2) override;
+ virtual void initGradients() override;
+ virtual void update() override;
virtual Cache &getVCache() const override;
- virtual Mem &getVMem() const override;
#endif
};
diff --git a/tbox/src/wHex8.inl b/tbox/src/wHex8.inl
index 6912e483..a2c43adf 100644
--- a/tbox/src/wHex8.inl
+++ b/tbox/src/wHex8.inl
@@ -14,14 +14,6 @@
* limitations under the License.
*/
-/**
- * @brief Return private memory for Hex8 private methods
- */
-inline MemHex8 &Hex8::getMem() const
-{
- return *pmem;
-}
-
/**
* @brief Return cache for Hex8 private methods
*/
diff --git a/tbox/src/wLine2.cpp b/tbox/src/wLine2.cpp
index 85dce33b..b28d1316 100644
--- a/tbox/src/wLine2.cpp
+++ b/tbox/src/wLine2.cpp
@@ -18,7 +18,6 @@
#include "wSfLine2.h"
#include "wGaussLine2.h"
#include "wCacheLine2.h"
-#include "wMemLine2.h"
#include "wNode.h"
using namespace tbox;
@@ -26,61 +25,114 @@ using namespace tbox;
Line2::Line2(int n, Tag *_ptag, Tag *_etag,
std::vector<Tag *> &_parts, std::vector<Node *> &nods) : Element(n, _ptag, _etag, _parts, nods)
{
- pmem = new MemLine2(*this, static_cast<int>(getCache(order).gauss.getN()));
}
-Line2::~Line2()
+
+/**
+ * @brief Return cache holding shape functions and Gauss points for all Line2 elements
+ */
+Cache &Line2::getVCache() const
{
- if (pmem)
- delete pmem;
+ return getCache(order);
}
/**
- * @brief Return private memory handling Jacobian for this Line2 element
+ * @brief Initialize precomputed values
*/
-Mem &Line2::getVMem() const
+void Line2::initValues(bool isvol, int ordr)
{
- return getMem();
+ // Set flags
+ isVol = isvol;
+ order = ordr;
+ hasVals = true;
+ // Update values
+ size_t ngp = getCache(order).getVGauss().getN();
+ Js.resize(ngp);
+ iJs.resize(ngp);
+ detJs.resize(ngp);
+ this->update();
}
-
/**
- * @brief Return cache holding shape functions and Gauss points for all Line2 elements
+ * @brief Initialize precomputed gradients
*/
-Cache &Line2::getVCache() const
+void Line2::initGradients()
{
- return getCache(order);
+ // Set flags
+ hasGrads = true;
+ // Update values
+ size_t ngp = getCache(order).getVGauss().getN();
+ gradJs.resize(ngp);
+ gradDetJs.resize(ngp);
+ this->update();
+}
+/**
+ * @brief Update precomputed values and gradients
+ */
+void Line2::update()
+{
+ // Update values
+ if (hasVals)
+ {
+ // Jacobian
+ if (isVol)
+ throw std::runtime_error("Line2::update(): the problem must always have a dimension less than Line2!\n");
+ // Line2 is never of same dimension as Problem (surface element)
+ this->buildSurfaceJ();
+ // Length and CG
+ this->computeV();
+ this->computeCg();
+ // Normal
+ this->computeNormal();
+ }
+ // Update gradients
+ if (hasGrads)
+ {
+ // Jacobian
+ if (isVol)
+ throw std::runtime_error("Line2::updateGradient(): the problem must always have a dimension less than Line2!\n");
+ // Line2 is never of same dimension as Problem (surface element)
+ this->buildGradientSurfaceJ();
+ // Length
+ this->buildGradientV(2);
+ // Normal
+ this->computeGradientNormal();
+ }
}
/**
- * @brief Compute Jacobian determinant at Gauss point k
+ * @brief Compute Jacobian determinant
*
* detJ = norm(x_i * dN_i)
*/
-double Line2::buildJ(size_t k) const
+void Line2::buildSurfaceJ()
{
- return this->buildTangents(k)[0].norm();
+ for (size_t k = 0; k < getCache(order).getVGauss().getN(); ++k)
+ detJs[k] = this->buildTangents(k)[0].norm();
}
/**
- * @brief Compute Jacobian determinant gradient with respect to each node coordinate, at Gauss point k
+ * @brief Compute Jacobian determinant gradient with respect to each node coordinate
*
* dDetJ[l = i * 2 + j] = detJ_l * dt0_l;
*/
-void Line2::buildGradientJ(size_t k, std::vector<double> &dDetJ) const
+void Line2::buildGradientSurfaceJ()
{
- Eigen::MatrixXd const &dsf = getCache(order).getDsf(k);
-
- // Tangent vectors and determinant
- Eigen::Vector3d detJ = this->buildTangents(k)[0].normalized();
-
- // Jacobian gradient
- dDetJ.resize(2 * 2); // 2 nodes and 2 dimensions (x, y)
- for (size_t i = 0; i < 2; ++i)
+ CacheLine2 &cache = getCache(order);
+ for (size_t k = 0; k < cache.getVGauss().getN(); ++k)
{
- for (size_t j = 0; j < 2; ++j)
+ // Shape functions, tangent vectors and determinant
+ Eigen::MatrixXd const &dsf = getCache(order).getDsf(k);
+ Eigen::Vector3d detJ = this->buildTangents(k)[0].normalized();
+
+ // Jacobian gradient
+ gradDetJs[k].resize(2 * 2); // 2 nodes and 2 dimensions (x, y)
+ for (size_t i = 0; i < 2; ++i)
{
- Eigen::Vector3d dt0 = Eigen::Vector3d::Zero();
- dt0(j) = dsf(0, i); // retain sf gradients associated to ith node for dimension j
- dDetJ[i * 2 + j] = detJ.dot(dt0);
+ for (size_t j = 0; j < 2; ++j)
+ {
+ Eigen::Vector3d dt0 = Eigen::Vector3d::Zero();
+ dt0(j) = dsf(0, i); // retain sf gradients associated to ith node for dimension j
+ gradDetJs[k][i * 2 + j] = detJ.dot(dt0);
+ }
}
}
}
@@ -88,38 +140,34 @@ void Line2::buildGradientJ(size_t k, std::vector<double> &dDetJ) const
/**
* @brief Compute the element length
*/
-double Line2::computeV() const
+void Line2::computeV()
{
CacheLine2 &cache = getCache(order);
GaussLine2 &gauss = cache.gauss;
- MemLine2 &mem = getMem();
// Gauss integration
- double V = 0.0;
+ vol = 0.0;
for (size_t k = 0; k < gauss.getN(); ++k)
- V += gauss.getW(k) * mem.getDetJ(k);
- return V;
+ vol += gauss.getW(k) * this->getDetJ(k);
}
/**
* @brief Compute length gradient with respect to each node coordinate
*
* dV[i] = w_k * dDetJ[i]_k
*/
-std::vector<double> Line2::buildGradientV(int dim) const
+void Line2::buildGradientV(int dim)
{
GaussLine2 &gauss = getCache(order).gauss;
- MemLine2 &mem = getMem();
// Gradient of volume
- std::vector<double> dVol(2 * dim, 0.); // 2 nodes
+ gradVol.resize(2 * dim, 0.); // 2 nodes
for (size_t k = 0; k < gauss.getN(); ++k)
{
- double w = gauss.getW(k); // GP weight
- std::vector<double> dDetJ = mem.getGradDetJ(k); // gradient of Jacobian determinant
+ double w = gauss.getW(k); // GP weight
+ std::vector<double> dDetJ = this->getGradDetJ(k); // gradient of Jacobian determinant
for (size_t i = 0; i < dDetJ.size(); ++i)
- dVol[i] += w * dDetJ[i];
+ gradVol[i] += w * dDetJ[i];
}
- return dVol;
}
/**
@@ -149,16 +197,16 @@ std::vector<Eigen::Vector3d> Line2::buildTangents(size_t k) const
/**
* @brief Compute element unit normal vector
*/
-Eigen::Vector3d Line2::normal() const
+void Line2::computeNormal()
{
Eigen::Vector3d t = this->computeTangents()[0].normalized();
- return Eigen::Vector3d(-t(1), t(0), 0.);
+ normal = Eigen::Vector3d(-t(1), t(0), 0.);
}
/**
* @brief Compute gradient of element unit normal vector with respect to each node coordinate
*/
-std::vector<Eigen::Vector3d> Line2::gradientNormal() const
+void Line2::computeGradientNormal()
{
// Compute normal and norm
Eigen::Vector3d t = this->computeTangents()[0];
@@ -167,10 +215,9 @@ std::vector<Eigen::Vector3d> Line2::gradientNormal() const
// Compute gradients
Eigen::Matrix3d inn = (Eigen::Matrix3d::Identity() - n * n.transpose() / (nn * nn)) / nn; // gradient of n/|n|
- std::vector<Eigen::Vector3d> dn(2 * 2); // 2 nodes and 2 dimensions
- dn[0] = inn * Eigen::Vector3d(0, -1, 0);
- dn[1] = inn * Eigen::Vector3d(1, 0, 0);
- dn[2] = inn * Eigen::Vector3d(0, 1, 0);
- dn[3] = inn * Eigen::Vector3d(-1, 0, 0); // outer gradient times gradient of t (pos) wtr x
- return dn;
+ gradNormal.resize(2 * 2); // 2 nodes and 2 dimensions
+ gradNormal[0] = inn * Eigen::Vector3d(0, -1, 0);
+ gradNormal[1] = inn * Eigen::Vector3d(1, 0, 0);
+ gradNormal[2] = inn * Eigen::Vector3d(0, 1, 0);
+ gradNormal[3] = inn * Eigen::Vector3d(-1, 0, 0); // outer gradient times gradient of t (pos) wtr x
}
diff --git a/tbox/src/wLine2.h b/tbox/src/wLine2.h
index 9e6e5f0d..6b076608 100644
--- a/tbox/src/wLine2.h
+++ b/tbox/src/wLine2.h
@@ -32,33 +32,29 @@ namespace tbox
class TBOX_API Line2 : public Element
{
#ifndef SWIG
- friend class MemLine2;
-
private:
- MemLine2 *pmem; ///< private memory of precalculated values
- inline MemLine2 &getMem() const;
inline static CacheLine2 &getCache(int n);
protected:
- virtual double buildJ(size_t k) const override;
- virtual void buildGradientJ(size_t k, std::vector<double> &dDetJ) const override;
-
-public: // @todo
- virtual double computeV() const override;
- virtual std::vector<double> buildGradientV(int dim) const override;
- virtual Eigen::Vector3d normal() const override;
- virtual std::vector<Eigen::Vector3d> gradientNormal() const override;
+ virtual void buildSurfaceJ() override;
+ virtual void buildGradientSurfaceJ() override;
+ virtual void computeV() override;
+ virtual void buildGradientV(int dim) override;
+ virtual void computeNormal() override;
+ virtual void computeGradientNormal() override;
virtual std::vector<Eigen::Vector3d> computeTangents() const override;
virtual std::vector<Eigen::Vector3d> buildTangents(size_t k) const override;
#endif
public:
Line2(int n, Tag *_ptag, Tag *_etag, std::vector<Tag *> &_parts, std::vector<Node *> &nods);
- virtual ~Line2();
+ virtual ~Line2() {}
virtual ELTYPE type() const override { return ELTYPE::LINE2; }
#ifndef SWIG
- virtual Mem &getVMem() const override;
+ virtual void initValues(bool isvol, int ordr = 2) override;
+ virtual void initGradients() override;
+ virtual void update() override;
virtual Cache &getVCache() const override;
#endif
};
diff --git a/tbox/src/wLine2.inl b/tbox/src/wLine2.inl
index 704f2314..5dc09184 100644
--- a/tbox/src/wLine2.inl
+++ b/tbox/src/wLine2.inl
@@ -14,14 +14,6 @@
* limitations under the License.
*/
-/**
- * @brief Return private memory for Line2 private methods
- */
-inline MemLine2 &Line2::getMem() const
-{
- return *pmem;
-}
-
/**
* @brief Return cache for Line2 private methods
*/
diff --git a/tbox/src/wMem.cpp b/tbox/src/wMem.cpp
deleted file mode 100644
index 4a003a0c..00000000
--- a/tbox/src/wMem.cpp
+++ /dev/null
@@ -1,59 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#include "wMem.h"
-#include "wElement.h"
-
-using namespace tbox;
-
-Mem::Mem(size_t _ngp) : valid(false), validG(false), ngp(_ngp), detJs(_ngp), Js(_ngp), Jinvs(_ngp)
-{
- vol = 0;
- cg = Eigen::Vector3d::Zero();
-}
-
-/**
- * @brief Check that determinant of Jacobian is positive
- */
-void Mem::checkJac(Element const &e, double detJ) const
-{
- if (detJ <= 0)
- {
- std::stringstream err;
- err << "Negative jacobian detected for element: " << e << "\n";
- throw std::runtime_error(err.str());
- }
-}
-
-/**
- * @brief Initialize the gradients
- */
-void Mem::initGradients(bool sameDim)
-{
- gradJs.resize(ngp);
- gradDetJs.resize(ngp);
- this->updateGradients(sameDim);
-}
-
-void Mem::update(bool sameDim)
-{
- throw std::runtime_error("Mem::update() not implemented!\n");
-}
-
-void Mem::updateGradients(bool sameDim)
-{
- throw std::runtime_error("Mem::updateGradients() not implemented!\n");
-}
diff --git a/tbox/src/wMem.h b/tbox/src/wMem.h
deleted file mode 100644
index 8a8945f4..00000000
--- a/tbox/src/wMem.h
+++ /dev/null
@@ -1,86 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#ifndef WMEM_H
-#define WMEM_H
-
-#include "tbox.h"
-#include <vector>
-#include <Eigen/Dense>
-
-#ifndef SWIG
-
-namespace tbox
-{
-
-/**
- * @brief Base class for element private memory
- * @authors Romain Boman, Adrien Crovato
- * @todo Remove capability to getJ and getJinv when not volume element
- * @todo Split values and gradients in 2 different Mem?
- */
-class TBOX_API Mem
-{
-protected:
- // Flag
- bool valid; ///< validity flag (true if memory has been updated once)
- bool validG; ///< validity flag for gradients (true if gradients have been initialized)
-
- // Element value
- size_t ngp; ///< number of Gauss points
- double vol; ///< element volume
- Eigen::Vector3d cg; ///< element center of gravity
-
- // GP values
- std::vector<double> detJs; ///< determinant of Jacobian matrix
- std::vector<Eigen::MatrixXd> Js; ///< Jacobian matrix
- std::vector<Eigen::MatrixXd> Jinvs; ///< inverse of Jacobian matrix
-
- // Gradients
- std::vector<double> gradVol; ///< gradient of volume wrt to nodes position
- std::vector<std::vector<double>> gradDetJs; ///< gradient of Jacobian determinant wrt to nodes position
- std::vector<std::vector<Eigen::MatrixXd>> gradJs; ///< gradient of Jacobian matrix wrt to nodes position
-
- // Functions
- void checkJac(Element const &e, double detJ) const;
- virtual void updateGradients(bool sameDim);
-
-public:
- explicit Mem(size_t _ngp);
- virtual ~Mem() {}
-
- inline bool isValid();
- inline bool isGradientValid();
- inline double getDetJ(size_t k);
- inline Eigen::MatrixXd const &getJ(size_t k);
- inline Eigen::MatrixXd const &getJinv(size_t k);
- inline double getVol();
- inline Eigen::Vector3d const &getCg();
- inline std::vector<double> const &getGradDetJ(size_t k);
- inline std::vector<Eigen::MatrixXd> const &getGradJ(size_t k);
- inline std::vector<double> const &getGradVol();
-
- void initGradients(bool sameDim);
- virtual void update(bool sameDim);
-};
-
-#include "wMem.inl"
-
-} // namespace tbox
-
-#endif
-
-#endif //WMEM_H
diff --git a/tbox/src/wMem.inl b/tbox/src/wMem.inl
deleted file mode 100644
index 217af82a..00000000
--- a/tbox/src/wMem.inl
+++ /dev/null
@@ -1,119 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-/**
- * @brief Return the validity flag
- */
-inline bool Mem::isValid()
-{
- return valid;
-}
-
-/**
- * @brief Return the validity flag (gradients)
- */
-inline bool Mem::isGradientValid()
-{
- return validG;
-}
-
-/**
- * @brief Return the determinant of the Jacobian at Gauss point k
- */
-inline double Mem::getDetJ(size_t k)
-{
- if (valid)
- return detJs[k];
- else
- throw std::runtime_error("Mem::getDetJ is invalid because element memory is not up-to-date!\n");
-}
-
-/**
- * @brief Return the the Jacobian matrix at Gauss point k
- */
-inline Eigen::MatrixXd const &Mem::getJ(size_t k)
-{
- if (valid)
- return Js[k];
- else
- throw std::runtime_error("Mem::getJ is invalid because element memory is not up-to-date!\n");
-}
-
-/**
- * @brief Return the inverse of the Jacobian matrix at Gauss point k
- */
-inline Eigen::MatrixXd const &Mem::getJinv(size_t k)
-{
- if (valid)
- return Jinvs[k];
- else
- throw std::runtime_error("Mem::getJinv is invalid because element memory is not up-to-date!\n");
-}
-
-/**
- * @brief Return the volume of the element
- */
-inline double Mem::getVol()
-{
- if (valid)
- return vol;
- else
- throw std::runtime_error("Mem::getVol is invalid because element memory is not up-to-date!\n");
-}
-
-/**
- * @brief Return the center of gravity of the element
- */
-inline Eigen::Vector3d const &Mem::getCg()
-{
- if (valid)
- return cg;
- else
- throw std::runtime_error("Mem::getCg is invalid because element memory is not up-to-date!\n");
-}
-
-/**
- * @brief Return the gradient of determinant of the Jacobian with respect to node positions at Gauss point k
- */
-inline std::vector<double> const &Mem::getGradDetJ(size_t k)
-{
- if (validG)
- return gradDetJs[k];
- else
- throw std::runtime_error("Mem::getGradDetJ is invalid because element memory (gradient) is not up-to-date!\n");
-}
-
-/**
- * @brief Return the gradient of the Jacobian matrix with respect to node positions at Gauss point k
- */
-inline std::vector<Eigen::MatrixXd> const &Mem::getGradJ(size_t k)
-{
- if (validG)
- return gradJs[k];
- else
- throw std::runtime_error("Mem::getGradJ is invalid because element memory (gradient) is not up-to-date!\n");
-}
-
-/**
- * @brief Return the gradient of the volume with respect to node positions
- */
-inline std::vector<double> const &Mem::getGradVol()
-{
- if (validG)
- return gradVol;
- else
- throw std::runtime_error("Mem::getGradVol is invalid because element memory (gradient) is not up-to-date!\n");
-}
diff --git a/tbox/src/wMemHex8.cpp b/tbox/src/wMemHex8.cpp
deleted file mode 100644
index c815f202..00000000
--- a/tbox/src/wMemHex8.cpp
+++ /dev/null
@@ -1,68 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#include "wMemHex8.h"
-#include "wNode.h"
-
-using namespace tbox;
-
-MemHex8::MemHex8(Hex8 &_e, size_t _ngp) : Mem(_ngp), e(_e)
-{
-}
-
-/**
- * @brief Update Hex8 values
- */
-void MemHex8::update(bool sameDim)
-{
- // Jacobian
- if (!sameDim)
- throw std::runtime_error("MemHex8::update(): the problem must always have the same dimension as Hex8!\n");
- // Hex8 is always of same dimension as Problem (volume element)
- for (int k = 0; k < ngp; ++k)
- {
- detJs[k] = e.buildJ(k, Js[k], Jinvs[k]);
- checkJac(e, detJs[k]);
- }
- // CG
- cg = Eigen::Vector3d::Zero();
- for (auto n : e.nodes)
- cg += n->pos;
- cg /= e.nodes.size();
- // Memory is now valid and the volume can be computed
- valid = true;
- vol = e.computeV();
-
- // Update gradient if they are needed
- if (validG)
- this->updateGradients(sameDim);
-}
-
-/**
- * @brief Update Hex8 gradients
- */
-void MemHex8::updateGradients(bool sameDim)
-{
- // Jacobian
- if (!sameDim)
- throw std::runtime_error("MemHex8::updateGradients(): the problem must always have the same dimension as Hex8!\n");
- // Hex8 is always of same dimension as Problem (volume element)
- for (int k = 0; k < ngp; ++k)
- e.buildGradientJ(k, gradDetJs[k], gradJs[k]);
- // Volume
- validG = true;
- gradVol = e.buildGradientV(3);
-}
diff --git a/tbox/src/wMemHex8.h b/tbox/src/wMemHex8.h
deleted file mode 100644
index 0787c7bb..00000000
--- a/tbox/src/wMemHex8.h
+++ /dev/null
@@ -1,50 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#ifndef WMEMHEX8_H
-#define WMEMHEX8_H
-
-#include "tbox.h"
-#include "wMem.h"
-#include "wHex8.h"
-
-#ifndef SWIG
-
-namespace tbox
-{
-
-/**
- * @brief Precalculated values for Hex8
- * @authors Romain Boman, Adrien Crovato
- */
-class TBOX_API MemHex8 : public Mem
-{
- // Element type
- Hex8 &e;
-
-public:
- explicit MemHex8(Hex8 &_e, size_t _ngp);
- virtual ~MemHex8() {}
-
- virtual void update(bool sameDim) override;
- virtual void updateGradients(bool sameDim) override;
-};
-
-} // namespace tbox
-
-#endif
-
-#endif //WMEMHEX8_H
diff --git a/tbox/src/wMemLine2.cpp b/tbox/src/wMemLine2.cpp
deleted file mode 100644
index 8e0c8022..00000000
--- a/tbox/src/wMemLine2.cpp
+++ /dev/null
@@ -1,65 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#include "wMemLine2.h"
-#include "wNode.h"
-
-using namespace tbox;
-
-MemLine2::MemLine2(Line2 &_e, size_t _ngp) : Mem(_ngp), e(_e)
-{
-}
-
-/**
- * @brief Update Line2 values
- */
-void MemLine2::update(bool sameDim)
-{
- // Jacobian
- if (sameDim)
- throw std::runtime_error("MemLine2::update(): the problem must always have a dimension less than Line2!\n");
- // Line2 is never of same dimension as Problem (surface element)
- for (int k = 0; k < ngp; ++k)
- detJs[k] = e.buildJ(k);
- // CG
- cg = Eigen::Vector3d::Zero();
- for (auto n : e.nodes)
- cg += n->pos;
- cg /= e.nodes.size();
- // Memory is now valid and the volume can be computed
- valid = true;
- vol = e.computeV();
-
- // Update gradient if they are needed
- if (validG)
- this->updateGradients(sameDim);
-}
-
-/**
- * @brief Update Line2 gradients
- */
-void MemLine2::updateGradients(bool sameDim)
-{
- // Jacobian
- if (sameDim)
- throw std::runtime_error("MemLine2::updateGradient(): the problem must always have a dimension less than Line2!\n");
- // Line2 is never of same dimension as Problem (surface element)
- for (int k = 0; k < ngp; ++k)
- e.buildGradientJ(k, gradDetJs[k]);
- // Volume
- validG = true;
- gradVol = e.buildGradientV(2);
-}
diff --git a/tbox/src/wMemLine2.h b/tbox/src/wMemLine2.h
deleted file mode 100644
index 38624af0..00000000
--- a/tbox/src/wMemLine2.h
+++ /dev/null
@@ -1,50 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#ifndef WMEMLINE2_H
-#define WMEMLINE2_H
-
-#include "tbox.h"
-#include "wMem.h"
-#include "wLine2.h"
-
-#ifndef SWIG
-
-namespace tbox
-{
-
-/**
- * @brief Precalculated values for Line2
- * @authors Romain Boman, Adrien Crovato
- */
-class TBOX_API MemLine2 : public Mem
-{
- // Element type
- Line2 &e;
-
-public:
- explicit MemLine2(Line2 &_e, size_t _ngp);
- virtual ~MemLine2() {}
-
- virtual void update(bool sameDim) override;
- virtual void updateGradients(bool sameDim) override;
-};
-
-} // namespace tbox
-
-#endif
-
-#endif //WMEMLINE2_H
diff --git a/tbox/src/wMemQuad4.cpp b/tbox/src/wMemQuad4.cpp
deleted file mode 100644
index 50ee3b51..00000000
--- a/tbox/src/wMemQuad4.cpp
+++ /dev/null
@@ -1,81 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#include "wMemQuad4.h"
-#include "wNode.h"
-
-using namespace tbox;
-
-MemQuad4::MemQuad4(Quad4 &_e, size_t _ngp) : Mem(_ngp), e(_e)
-{
-}
-
-/**
- * @brief Update Quad4 values
- */
-void MemQuad4::update(bool sameDim)
-{
- // Jacobian
- // Quad4 is not of same dimension as Problem (surface element)
- if (!sameDim)
- {
- for (int k = 0; k < ngp; ++k)
- detJs[k] = e.buildJ(k);
- }
- // Quad4 is of same dimension as Problem (volume element)
- else
- {
- for (int k = 0; k < ngp; ++k)
- {
- detJs[k] = e.buildJ(k, Js[k], Jinvs[k]);
- checkJac(e, detJs[k]);
- }
- }
- // CG
- cg = Eigen::Vector3d::Zero();
- for (auto n : e.nodes)
- cg += n->pos;
- cg /= e.nodes.size();
- // Memory is now valid and the volume can be computed
- valid = true;
- vol = e.computeV();
-
- // Update gradient if they are needed
- if (validG)
- this->updateGradients(sameDim);
-}
-
-/**
- * @brief Update Quad4 gradients
- */
-void MemQuad4::updateGradients(bool sameDim)
-{
- // Jacobian
- // Quad4 is not of same dimension as Problem (surface element)
- if (!sameDim)
- for (int k = 0; k < ngp; ++k)
- e.buildGradientJ(k, gradDetJs[k]);
- // Quad4 is of same dimension as Problem (volume element)
- else
- for (int k = 0; k < ngp; ++k)
- e.buildGradientJ(k, gradDetJs[k], gradJs[k]);
- // Volume
- validG = true;
- if (!sameDim)
- gradVol = e.buildGradientV(3); // gradient of surface
- else
- gradVol = e.buildGradientV(2); // gradient of volume
-}
diff --git a/tbox/src/wMemQuad4.h b/tbox/src/wMemQuad4.h
deleted file mode 100644
index cfe653f1..00000000
--- a/tbox/src/wMemQuad4.h
+++ /dev/null
@@ -1,50 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#ifndef WMEMQUAD4_H
-#define WMEMQUAD4_H
-
-#include "tbox.h"
-#include "wMem.h"
-#include "wQuad4.h"
-
-#ifndef SWIG
-
-namespace tbox
-{
-
-/**
- * @brief Precalculated values for Quad4
- * @authors Romain Boman, Adrien Crovato
- */
-class TBOX_API MemQuad4 : public Mem
-{
- // Element type
- Quad4 &e;
-
-public:
- explicit MemQuad4(Quad4 &_e, size_t _ngp);
- virtual ~MemQuad4() {}
-
- virtual void update(bool sameDim) override;
- virtual void updateGradients(bool sameDim) override;
-};
-
-} // namespace tbox
-
-#endif
-
-#endif //WMEMQUAD4_H
diff --git a/tbox/src/wMemTetra4.cpp b/tbox/src/wMemTetra4.cpp
deleted file mode 100644
index 5be98b0c..00000000
--- a/tbox/src/wMemTetra4.cpp
+++ /dev/null
@@ -1,68 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#include "wMemTetra4.h"
-#include "wNode.h"
-
-using namespace tbox;
-
-MemTetra4::MemTetra4(Tetra4 &_e, size_t _ngp) : Mem(_ngp), e(_e)
-{
-}
-
-/**
- * @brief Update Tetra4 values
- */
-void MemTetra4::update(bool sameDim)
-{
- // Jacobian
- if (!sameDim)
- throw std::runtime_error("MemTetra4::update(): the problem must always have the same dimension as Tetra4!\n");
- // Tetra4 is always of same dimension as Problem (volume element)
- for (size_t k = 0; k < ngp; ++k)
- {
- detJs[k] = e.buildJ(k, Js[k], Jinvs[k]);
- checkJac(e, detJs[k]);
- }
- // CG
- cg = Eigen::Vector3d::Zero();
- for (auto n : e.nodes)
- cg += n->pos;
- cg /= e.nodes.size();
- // Memory is now valid and the volume can be computed
- valid = true;
- vol = e.computeV();
-
- // Update gradient if they are needed
- if (validG)
- this->updateGradients(sameDim);
-}
-
-/**
- * @brief Update Tetra4 gradients
- */
-void MemTetra4::updateGradients(bool sameDim)
-{
- // Jacobian
- if (!sameDim)
- throw std::runtime_error("MemTetra4::updateGradients(): the problem must always have the same dimension as Tetra4!\n");
- // Tetra4 is always of same dimension as Problem (volume element)
- for (int k = 0; k < ngp; ++k)
- e.buildGradientJ(k, gradDetJs[k], gradJs[k]);
- // Volume
- validG = true;
- gradVol = e.buildGradientV(3);
-}
diff --git a/tbox/src/wMemTetra4.h b/tbox/src/wMemTetra4.h
deleted file mode 100644
index 42a71dcb..00000000
--- a/tbox/src/wMemTetra4.h
+++ /dev/null
@@ -1,50 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#ifndef WMEMTETRA4_H
-#define WMEMTETRA4_H
-
-#include "tbox.h"
-#include "wMem.h"
-#include "wTetra4.h"
-
-#ifndef SWIG
-
-namespace tbox
-{
-
-/**
- * @brief Precalculated values for Tetra4
- * @authors Romain Boman, Adrien Crovato
- */
-class TBOX_API MemTetra4 : public Mem
-{
- // Element type
- Tetra4 &e;
-
-public:
- explicit MemTetra4(Tetra4 &_e, size_t _ngp);
- virtual ~MemTetra4() {}
-
- virtual void update(bool sameDim) override;
- virtual void updateGradients(bool sameDim) override;
-};
-
-} // namespace tbox
-
-#endif
-
-#endif //WMEMTETRA4_H
diff --git a/tbox/src/wMemTri3.cpp b/tbox/src/wMemTri3.cpp
deleted file mode 100644
index 707e7d55..00000000
--- a/tbox/src/wMemTri3.cpp
+++ /dev/null
@@ -1,81 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#include "wMemTri3.h"
-#include "wNode.h"
-
-using namespace tbox;
-
-MemTri3::MemTri3(Tri3 &_e, size_t _ngp) : Mem(_ngp), e(_e)
-{
-}
-
-/**
- * @brief Update Tri3 values
- */
-void MemTri3::update(bool sameDim)
-{
- // Jacobian
- // Tri3 is not of same dimension as Problem (surface element)
- if (!sameDim)
- {
- for (int k = 0; k < ngp; ++k)
- detJs[k] = e.buildJ(k);
- }
- // Tri3 is of same dimension as Problem (volume element)
- else
- {
- for (int k = 0; k < ngp; ++k)
- {
- detJs[k] = e.buildJ(k, Js[k], Jinvs[k]);
- checkJac(e, detJs[k]);
- }
- }
- // CG
- cg = Eigen::Vector3d::Zero();
- for (auto n : e.nodes)
- cg += n->pos;
- cg /= e.nodes.size();
- // Memory is now valid and the volume can be computed
- valid = true;
- vol = e.computeV();
-
- // Update gradient if they are needed
- if (validG)
- this->updateGradients(sameDim);
-}
-
-/**
- * @brief Update Tri3 gradients
- */
-void MemTri3::updateGradients(bool sameDim)
-{
- // Jacobian
- // Tri3 is not of same dimension as Problem (surface element)
- if (!sameDim)
- for (int k = 0; k < ngp; ++k)
- e.buildGradientJ(k, gradDetJs[k]);
- // Tri3 is of same dimension as Problem (volume element)
- else
- for (int k = 0; k < ngp; ++k)
- e.buildGradientJ(k, gradDetJs[k], gradJs[k]);
- // Volume
- validG = true;
- if (!sameDim)
- gradVol = e.buildGradientV(3); // gradient of surface
- else
- gradVol = e.buildGradientV(2); // gradient of volume
-}
diff --git a/tbox/src/wMemTri3.h b/tbox/src/wMemTri3.h
deleted file mode 100644
index 025e54ea..00000000
--- a/tbox/src/wMemTri3.h
+++ /dev/null
@@ -1,50 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#ifndef WMEMTRI3_H
-#define WMEMTRI3_H
-
-#include "tbox.h"
-#include "wMem.h"
-#include "wTri3.h"
-
-#ifndef SWIG
-
-namespace tbox
-{
-
-/**
- * @brief Precalculated values for Tri3
- * @authors Romain Boman, Adrien Crovato
- */
-class TBOX_API MemTri3 : public Mem
-{
- // Element type
- Tri3 &e;
-
-public:
- explicit MemTri3(Tri3 &_e, size_t _ngp);
- virtual ~MemTri3() {}
-
- virtual void update(bool sameDim) override;
- virtual void updateGradients(bool sameDim) override;
-};
-
-} // namespace tbox
-
-#endif
-
-#endif //WMEMTRI3_H
diff --git a/tbox/src/wMshCrack.cpp b/tbox/src/wMshCrack.cpp
index a54dc692..8b550193 100644
--- a/tbox/src/wMshCrack.cpp
+++ b/tbox/src/wMshCrack.cpp
@@ -225,11 +225,10 @@ void MshCrack::modifyBoundaries()
std::vector<Eigen::Vector3d> crkCG(crckGrp->tag->elems.size()), crkN(crckGrp->tag->elems.size());
for (size_t i = 0; i < crckGrp->tag->elems.size(); ++i)
{
- crkCG[i] = Eigen::Vector3d::Zero();
- for (size_t j = 0; j < crckGrp->tag->elems[i]->nodes.size(); ++j)
- crkCG[i] += crckGrp->tag->elems[i]->nodes[j]->pos;
- crkCG[i] /= crckGrp->tag->elems[i]->nodes.size();
- crkN[i] = crckGrp->tag->elems[i]->normal();
+ crckGrp->tag->elems[i]->computeCg();
+ crkCG[i] = crckGrp->tag->elems[i]->cg;
+ crckGrp->tag->elems[i]->computeNormal();
+ crkN[i] = crckGrp->tag->elems[i]->normal;
}
// find elements on the opposite side of the crack and swap their nodes
for (auto g : bndGrp)
diff --git a/tbox/src/wMshDeform.cpp b/tbox/src/wMshDeform.cpp
index 6eacfeba..259fa279 100644
--- a/tbox/src/wMshDeform.cpp
+++ b/tbox/src/wMshDeform.cpp
@@ -23,7 +23,6 @@
#include "wElement.h"
#include "wCache.h"
#include "wGauss.h"
-#include "wMem.h"
#include "wGmres.h"
#include <deque>
@@ -241,7 +240,6 @@ Eigen::MatrixXd MshDeform::buildK(tbox::Element const &e, Eigen::MatrixXd const
// Get precomputed values
Cache &cache = e.getVCache();
Gauss &gauss = cache.getVGauss();
- Mem &mem = e.getVMem();
size_t ns = e.nodes.size(); // number of nodes
// Elementary stiffness matrix
@@ -249,7 +247,7 @@ Eigen::MatrixXd MshDeform::buildK(tbox::Element const &e, Eigen::MatrixXd const
for (size_t k = 0; k < gauss.getN(); ++k)
{
// Jacobian inverse
- Eigen::MatrixXd const &iJ = mem.getJinv(k);
+ Eigen::MatrixXd const &iJ = e.getJinv(k);
// Fill B matrix (shape functions)
Eigen::MatrixXd const &dff = cache.getDsf(k);
Eigen::MatrixXd B;
@@ -286,7 +284,7 @@ Eigen::MatrixXd MshDeform::buildK(tbox::Element const &e, Eigen::MatrixXd const
}
// Compute stiffness matrix
- K += B.transpose() * H * B * gauss.getW(k) * mem.getDetJ(k);
+ K += B.transpose() * H * B * gauss.getW(k) * e.getDetJ(k);
}
return K;
}
@@ -314,8 +312,8 @@ void MshDeform::build(Eigen::SparseMatrix<double, Eigen::RowMajor> &K)
if (nDim == 2)
{
H = Eigen::Matrix3d::Zero();
- double nu = 0.; // Poisson ratio [-1, 0.5]
- double E = 1 / e->computeV(); // Young modulus
+ double nu = 0.; // Poisson ratio [-1, 0.5]
+ double E = 1 / e->getVol(); // Young modulus
H(0, 0) = 1 - nu;
H(1, 1) = 1 - nu;
H(2, 2) = (1 - 2 * nu) / 2;
@@ -326,8 +324,8 @@ void MshDeform::build(Eigen::SparseMatrix<double, Eigen::RowMajor> &K)
else
{
H = Eigen::MatrixXd::Zero(6, 6);
- double nu = 0.; // Poisson ratio [-1, 0.5]
- double E = 1 / e->computeV(); // Young modulus
+ double nu = 0.; // Poisson ratio [-1, 0.5]
+ double E = 1 / e->getVol(); // Young modulus
H(0, 0) = 1 - nu;
H(1, 1) = 1 - nu;
H(2, 2) = 1 - nu;
@@ -494,13 +492,13 @@ void MshDeform::deform()
}
// Update element memory
for (auto e : fldElems)
- e->getVMem().update(true);
+ e->update();
for (auto e : movBndElems)
- e->getVMem().update(false);
+ e->update();
for (auto e : symBndElems)
- e->getVMem().update(false);
+ e->update();
for (auto e : intBndElems)
- e->getVMem().update(false);
+ e->update();
std::cout << "done" << std::endl;
}
diff --git a/tbox/src/wQuad4.cpp b/tbox/src/wQuad4.cpp
index d3981e8b..c2ee84f6 100644
--- a/tbox/src/wQuad4.cpp
+++ b/tbox/src/wQuad4.cpp
@@ -18,7 +18,6 @@
#include "wSfQuad4.h"
#include "wGaussQuad4.h"
#include "wCacheQuad4.h"
-#include "wMemQuad4.h"
#include "wNode.h"
using namespace tbox;
@@ -26,72 +25,135 @@ using namespace tbox;
Quad4::Quad4(int n, Tag *_ptag, Tag *_etag,
std::vector<Tag *> &_parts, std::vector<Node *> &nods) : Element(n, _ptag, _etag, _parts, nods)
{
- pmem = new MemQuad4(*this, static_cast<int>(getCache(order).gauss.getN()));
}
-Quad4::~Quad4()
+
+/**
+ * @brief Return cache holding shape functions and Gauss points for all Quad4 elements
+ */
+Cache &Quad4::getVCache() const
{
- if (pmem)
- delete pmem;
+ return getCache(order);
}
/**
- * @brief Return private memory handling Jacobian for this Quad4 element
+ * @brief Initialize precomputed values
*/
-Mem &Quad4::getVMem() const
+void Quad4::initValues(bool isvol, int ordr)
{
- return getMem();
+ // Set flags
+ isVol = isvol;
+ order = ordr;
+ hasVals = true;
+ // Update values
+ size_t ngp = getCache(order).getVGauss().getN();
+ Js.resize(ngp);
+ iJs.resize(ngp);
+ detJs.resize(ngp);
+ this->update();
}
-
/**
- * @brief Return cache holding shape functions and Gauss points for all Quad4 elements
+ * @brief Initialize precomputed gradients
*/
-Cache &Quad4::getVCache() const
+void Quad4::initGradients()
{
- return getCache(order);
+ // Set flags
+ hasGrads = true;
+ // Update values
+ size_t ngp = getCache(order).getVGauss().getN();
+ gradJs.resize(ngp);
+ gradDetJs.resize(ngp);
+ this->update();
+}
+/**
+ * @brief Update precomputed values and gradients
+ */
+void Quad4::update()
+{
+ // Update values
+ if (hasVals)
+ {
+ // Jacobian and normal
+ // Quad4 is not of same dimension as Problem (surface element)
+ if (!isVol)
+ {
+ this->buildSurfaceJ();
+ this->computeNormal();
+ }
+ // Quad4 is of same dimension as Problem (volume element)
+ else
+ this->buildJ();
+ // Volume and CG
+ this->computeV();
+ this->computeCg();
+ }
+ // Update gradients
+ if (hasGrads)
+ {
+ // Jacobian, volume and normal
+ // Quad4 is not of same dimension as Problem (surface element)
+ if (!isVol)
+ {
+ this->buildGradientSurfaceJ();
+ this->buildGradientV(3); // gradient of surface
+ this->computeGradientNormal();
+ }
+ // Quad4 is of same dimension as Problem (volume element)
+ else
+ {
+ this->buildGradientJ();
+ this->buildGradientV(2); // gradient of volume
+ }
+ }
}
/**
- * @brief Compute Jacobian determinant at Gauss point k
+ * @brief Compute Jacobian determinant
*
* detJ = norm(x_i dN_i,xi X x_i dN_i,eta)
*/
-double Quad4::buildJ(size_t k) const
+void Quad4::buildSurfaceJ()
{
- std::vector<Eigen::Vector3d> ts = this->buildTangents(k);
- return ts[0].cross(ts[1]).norm();
+ for (size_t k = 0; k < getCache(order).getVGauss().getN(); ++k)
+ {
+ std::vector<Eigen::Vector3d> ts = this->buildTangents(k);
+ detJs[k] = ts[0].cross(ts[1]).norm();
+ }
}
/**
* @brief Compute Jacobian matrix, inverse and determinant at Gauss point k
*
* J_ij(2,2) = diN_n * xj_n (n is the node)
*/
-double Quad4::buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const
+void Quad4::buildJ()
{
- Eigen::MatrixXd const &dff = getCache(order).getDsf(k);
+ CacheQuad4 &cache = getCache(order);
+ for (size_t k = 0; k < cache.getVGauss().getN(); ++k)
+ {
+ Eigen::MatrixXd const &dff = cache.getDsf(k);
- Eigen::Matrix2d JJ = Eigen::Matrix2d::Zero(); // temporary fixed-size matrix to efficiently compute the inverse
- size_t i = 0;
- for (auto it = nodes.begin(); it != nodes.end(); ++it, ++i)
- JJ += dff.col(i) * (*it)->pos.transpose().block<1, 2>(0, 0);
- J = JJ;
- iJ = JJ.inverse();
- return JJ.determinant();
+ Eigen::Matrix2d JJ = Eigen::Matrix2d::Zero(); // temporary fixed-size matrix to efficiently compute the inverse
+ size_t i = 0;
+ for (auto it = nodes.begin(); it != nodes.end(); ++it, ++i)
+ JJ += dff.col(i) * (*it)->pos.transpose().block<1, 2>(0, 0);
+ Js[k] = JJ;
+ iJs[k] = JJ.inverse();
+ detJs[k] = JJ.determinant();
+ this->checkJac(k);
+ }
}
/**
* @brief Return the element surface
*/
-double Quad4::computeV() const
+void Quad4::computeV()
{
CacheQuad4 &cache = getCache(order);
GaussQuad4 &gauss = cache.gauss;
- MemQuad4 &mem = getMem();
// Gauss integration
- double V = 0.0;
+ vol = 0.0;
for (size_t k = 0; k < gauss.getN(); ++k)
- V += gauss.getW(k) * mem.getDetJ(k);
- return V;
+ vol += gauss.getW(k) * this->getDetJ(k);
}
/**
@@ -126,8 +188,8 @@ std::vector<Eigen::Vector3d> Quad4::buildTangents(size_t k) const
/**
* @brief Compute element unit normal vector
*/
-Eigen::Vector3d Quad4::normal() const
+void Quad4::computeNormal()
{
std::vector<Eigen::Vector3d> ts = this->computeTangents();
- return ts[0].cross(ts[1]).normalized();
+ normal = ts[0].cross(ts[1]).normalized();
}
diff --git a/tbox/src/wQuad4.h b/tbox/src/wQuad4.h
index 0e3975b5..75cde412 100644
--- a/tbox/src/wQuad4.h
+++ b/tbox/src/wQuad4.h
@@ -32,32 +32,28 @@ namespace tbox
class TBOX_API Quad4 : public Element
{
#ifndef SWIG
- friend class MemQuad4;
-
private:
- MemQuad4 *pmem; ///< private memory of precalculated values
inline static CacheQuad4 &getCache(int n);
- inline MemQuad4 &getMem() const;
protected:
- virtual double buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const override;
- virtual double buildJ(size_t k) const override;
-
-public: // @todo
- virtual double computeV() const override;
+ virtual void buildJ() override;
+ virtual void buildSurfaceJ() override;
+ virtual void computeV() override;
virtual std::vector<Eigen::Vector3d> computeTangents() const override;
virtual std::vector<Eigen::Vector3d> buildTangents(size_t k) const override;
- virtual Eigen::Vector3d normal() const override;
+ virtual void computeNormal() override;
#endif
public:
Quad4(int n, Tag *_ptag, Tag *_etag, std::vector<Tag *> &_parts, std::vector<Node *> &nods);
- virtual ~Quad4();
+ virtual ~Quad4() {}
virtual ELTYPE type() const override { return ELTYPE::QUAD4; }
#ifndef SWIG
+ virtual void initValues(bool isvol, int ordr = 2) override;
+ virtual void initGradients() override;
+ virtual void update() override;
virtual Cache &getVCache() const override;
- virtual Mem &getVMem() const override;
#endif
};
diff --git a/tbox/src/wQuad4.inl b/tbox/src/wQuad4.inl
index 6e23c1b1..e77859ff 100644
--- a/tbox/src/wQuad4.inl
+++ b/tbox/src/wQuad4.inl
@@ -14,14 +14,6 @@
* limitations under the License.
*/
-/**
- * @brief Return private memory for Quad4 private methods
- */
-inline MemQuad4 &Quad4::getMem() const
-{
- return *pmem;
-}
-
/**
* @brief Return cache for Quad4 private methods
*/
diff --git a/tbox/src/wTetra4.cpp b/tbox/src/wTetra4.cpp
index e7b8e5e3..8c1b42ab 100644
--- a/tbox/src/wTetra4.cpp
+++ b/tbox/src/wTetra4.cpp
@@ -18,7 +18,6 @@
#include "wSfTetra4.h"
#include "wGaussTetra4.h"
#include "wCacheTetra4.h"
-#include "wMemTetra4.h"
#include "wNode.h"
using namespace tbox;
@@ -26,111 +25,160 @@ using namespace tbox;
Tetra4::Tetra4(int n, Tag *_ptag, Tag *_etag,
std::vector<Tag *> &_parts, std::vector<Node *> &nods) : Element(n, _ptag, _etag, _parts, nods)
{
- pmem = new MemTetra4(*this, static_cast<int>(getCache(order).gauss.getN()));
}
-Tetra4::~Tetra4()
+
+/**
+ * @brief Return cache handling shape functions and Gauss points for all Tetra4 elements
+ */
+Cache &Tetra4::getVCache() const
{
- if (pmem)
- delete pmem;
+ return getCache(order);
}
/**
- * @brief Return private memory handling Jacobian for this Tetra4 element
+ * @brief Initialize precomputed values
*/
-Mem &Tetra4::getVMem() const
+void Tetra4::initValues(bool isvol, int ordr)
{
- return getMem();
+ // Set flags
+ isVol = isvol;
+ order = ordr;
+ hasVals = true;
+ // Update values
+ size_t ngp = getCache(order).getVGauss().getN();
+ Js.resize(ngp);
+ iJs.resize(ngp);
+ detJs.resize(ngp);
+ this->update();
}
-
/**
- * @brief Return cache handling shape functions and Gauss points for all Tetra4 elements
+ * @brief Initialize precomputed gradients
*/
-Cache &Tetra4::getVCache() const
+void Tetra4::initGradients()
{
- return getCache(order);
+ // Set flags
+ hasGrads = true;
+ // Update values
+ size_t ngp = getCache(order).getVGauss().getN();
+ gradJs.resize(ngp);
+ gradDetJs.resize(ngp);
+ this->update();
+}
+/**
+ * @brief Update precomputed values and gradients
+ */
+void Tetra4::update()
+{
+ // Update gradients
+ if (hasVals)
+ {
+ // Jacobian
+ if (!isVol)
+ throw std::runtime_error("Tetra4::update(): the problem must always have the same dimension as Tetra4!\n");
+ // Tetra4 is always of same dimension as Problem (volume element)
+ this->buildJ();
+ // Volume and CG
+ this->computeV();
+ this->computeCg();
+ }
+ // Update gradients
+ if (hasGrads)
+ {
+ // Jacobian
+ if (!isVol)
+ throw std::runtime_error("Tetra4::updateGradients(): the problem must always have the same dimension as Tetra4!\n");
+ // Tetra4 is always of same dimension as Problem (volume element)
+ this->buildGradientJ();
+ // Volume
+ this->buildGradientV(3);
+ }
}
/**
- * @brief Compute Jacobian matrix, inverse and determinant at Gauss point k
+ * @brief Compute Jacobian matrix, inverse and determinant
*
* J_ij(3,3) = diN_n * xj_n (n is the node)
*/
-double Tetra4::buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const
+void Tetra4::buildJ()
{
- Eigen::MatrixXd const &dff = getCache(order).getDsf(k);
+ CacheTetra4 &cache = getCache(order);
+ for (size_t k = 0; k < cache.getVGauss().getN(); ++k)
+ {
+ Eigen::MatrixXd const &dff = cache.getDsf(k);
- Eigen::Matrix3d JJ = Eigen::Matrix3d::Zero(); // temporary fixed-size matrix to efficiently compute the inverse
- size_t i = 0;
- for (auto it = nodes.begin(); it != nodes.end(); ++it, ++i)
- JJ += dff.col(i) * ((*it)->pos).transpose();
- J = JJ;
- iJ = JJ.inverse();
- return JJ.determinant();
+ Eigen::Matrix3d JJ = Eigen::Matrix3d::Zero(); // temporary fixed-size matrix to efficiently compute the inverse
+ size_t i = 0;
+ for (auto it = nodes.begin(); it != nodes.end(); ++it, ++i)
+ JJ += dff.col(i) * ((*it)->pos).transpose();
+ Js[k] = JJ;
+ iJs[k] = JJ.inverse();
+ detJs[k] = JJ.determinant();
+ this->checkJac(k);
+ }
}
/**
- * @brief Compute Jacobian matrix gradient with respect to each node coordinate, at Gauss point k
+ * @brief Compute Jacobian matrix and determinant gradient with respect to each node coordinate
*
* dJ[l = i * 3 + j](3,3) = dN_i(3,1) * dxy_i(1,3)/dx_l (for node i and dimension j)
*/
-void Tetra4::buildGradientJ(size_t k, std::vector<double> &dDetJ, std::vector<Eigen::MatrixXd> &dJ) const
+void Tetra4::buildGradientJ()
{
- MemTetra4 &mem = getMem();
- Eigen::MatrixXd const &dsf = getCache(order).getDsf(k);
-
- // Jacobian
- dJ.resize(4 * 3); // 4 nodes and 3 dimensions (x, y, z)
- for (size_t i = 0; i < 4; ++i)
+ CacheTetra4 &cache = getCache(order);
+ for (size_t k = 0; k < cache.getVGauss().getN(); ++k)
{
- for (size_t j = 0; j < 3; ++j)
+ Eigen::MatrixXd const &dsf = getCache(order).getDsf(k);
+
+ // Jacobian
+ gradJs[k].resize(4 * 3); // 4 nodes and 3 dimensions (x, y, z)
+ for (size_t i = 0; i < 4; ++i)
{
- dJ[i * 3 + j] = Eigen::Matrix3d::Zero();
- dJ[i * 3 + j].col(j) = dsf.col(i); // retain sf gradients associated to ith node for dimension j
+ for (size_t j = 0; j < 3; ++j)
+ {
+ gradJs[k][i * 3 + j] = Eigen::Matrix3d::Zero();
+ gradJs[k][i * 3 + j].col(j) = dsf.col(i); // retain sf gradients associated to ith node for dimension j
+ }
}
+ // Determinant
+ double detJ = this->getDetJ(k);
+ Eigen::MatrixXd const &iJ = this->getJinv(k);
+ gradDetJs[k].resize(4 * 3);
+ for (size_t i = 0; i < 4 * 3; ++i)
+ gradDetJs[k][i] = detJ * (iJ * gradJs[k][i]).trace();
}
- // Determinant
- double detJ = mem.getDetJ(k);
- Eigen::MatrixXd const &iJ = mem.getJinv(k);
- dDetJ.resize(4 * 3);
- for (size_t i = 0; i < dDetJ.size(); ++i)
- dDetJ[i] = detJ * (iJ * dJ[i]).trace();
}
/**
* @brief Return the element volume
*/
-double Tetra4::computeV() const
+void Tetra4::computeV()
{
CacheTetra4 &cache = getCache(order);
GaussTetra4 &gauss = cache.gauss;
- MemTetra4 &mem = getMem();
// Gauss integration
- double V = 0.0;
+ vol = 0.0;
for (size_t k = 0; k < gauss.getN(); ++k)
- V += gauss.getW(k) * mem.getDetJ(k);
- return V;
+ vol += gauss.getW(k) * this->getDetJ(k);
}
/**
* @brief Compute volume gradient with respect to each node coordinate
*
* dV[i] = w_k * dDetJ[i]_k
*/
-std::vector<double> Tetra4::buildGradientV(int dim) const
+void Tetra4::buildGradientV(int dim)
{
GaussTetra4 &gauss = getCache(order).gauss;
- MemTetra4 &mem = getMem();
// Gradient of volume
- std::vector<double> dVol(4 * dim, 0.); // 4 nodes
+ gradVol.resize(4 * dim, 0.); // 4 nodes
for (size_t k = 0; k < gauss.getN(); ++k)
{
- double w = gauss.getW(k); // GP weight
- std::vector<double> dDetJ = mem.getGradDetJ(k); // gradient of Jacobian determinant
+ double w = gauss.getW(k); // GP weight
+ std::vector<double> dDetJ = this->getGradDetJ(k); // gradient of Jacobian determinant
for (size_t i = 0; i < dDetJ.size(); ++i)
- dVol[i] += w * dDetJ[i];
+ gradVol[i] += w * dDetJ[i];
}
- return dVol;
}
/**
@@ -143,5 +191,5 @@ Eigen::VectorXd Tetra4::computeGradient(std::vector<double> const &u, size_t k)
for (size_t i = 0; i < nodes.size(); ++i)
ue(i) = u[nodes[i]->row];
// gradient in global axis: inv(J)_ij * d_jN_i * ue_i
- return getMem().getJinv(k) * getCache(order).getDsf(k) * ue;
+ return this->getJinv(k) * getCache(order).getDsf(k) * ue;
}
diff --git a/tbox/src/wTetra4.h b/tbox/src/wTetra4.h
index cc2bdb39..371dcba8 100644
--- a/tbox/src/wTetra4.h
+++ b/tbox/src/wTetra4.h
@@ -32,30 +32,26 @@ namespace tbox
class TBOX_API Tetra4 : public Element
{
#ifndef SWIG
- friend class MemTetra4;
-
private:
- MemTetra4 *pmem; ///< private memory of precalculated values
inline static CacheTetra4 &getCache(int n);
- inline MemTetra4 &getMem() const;
protected:
- virtual double buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const override;
- virtual void buildGradientJ(size_t k, std::vector<double> &dDetJ, std::vector<Eigen::MatrixXd> &dJ) const override;
-
-public: // @todo
- virtual double computeV() const override;
- virtual std::vector<double> buildGradientV(int dim) const override;
+ virtual void buildJ() override;
+ virtual void buildGradientJ() override;
+ virtual void computeV() override;
+ virtual void buildGradientV(int dim) override;
#endif
public:
Tetra4(int n, Tag *_ptag, Tag *_etag, std::vector<Tag *> &_parts, std::vector<Node *> &nods);
- virtual ~Tetra4();
+ virtual ~Tetra4() {}
virtual ELTYPE type() const override { return ELTYPE::TETRA4; }
#ifndef SWIG
+ virtual void initValues(bool isvol, int ordr = 2) override;
+ virtual void initGradients() override;
+ virtual void update() override;
virtual Cache &getVCache() const override;
- virtual Mem &getVMem() const override;
virtual Eigen::VectorXd computeGradient(std::vector<double> const &u, size_t k) const override;
#endif
};
diff --git a/tbox/src/wTetra4.inl b/tbox/src/wTetra4.inl
index 506df073..64d5d52a 100644
--- a/tbox/src/wTetra4.inl
+++ b/tbox/src/wTetra4.inl
@@ -14,14 +14,6 @@
* limitations under the License.
*/
-/**
- * @brief Return private memory for Tetra4 private methods
- */
-inline MemTetra4 &Tetra4::getMem() const
-{
- return *pmem;
-}
-
/**
* @brief Return cache for Tetra4 private methods
*/
diff --git a/tbox/src/wTri3.cpp b/tbox/src/wTri3.cpp
index f145c4ed..ea16dd09 100644
--- a/tbox/src/wTri3.cpp
+++ b/tbox/src/wTri3.cpp
@@ -18,7 +18,6 @@
#include "wSfTri3.h"
#include "wGaussTri3.h"
#include "wCacheTri3.h"
-#include "wMemTri3.h"
#include "wNode.h"
using namespace tbox;
@@ -26,147 +25,214 @@ using namespace tbox;
Tri3::Tri3(int n, Tag *_ptag, Tag *_etag,
std::vector<Tag *> &_parts, std::vector<Node *> &nods) : Element(n, _ptag, _etag, _parts, nods)
{
- pmem = new MemTri3(*this, static_cast<int>(getCache(order).gauss.getN()));
}
-Tri3::~Tri3()
+
+/**
+ * @brief Return cache holding shape functions and Gauss points for all Tri3 elements
+ */
+Cache &Tri3::getVCache() const
{
- if (pmem)
- delete pmem;
+ return getCache(order);
}
/**
- * @brief Return private memory handling Jacobian for this Tri3 element
+ * @brief Initialize precomputed values
*/
-Mem &Tri3::getVMem() const
+void Tri3::initValues(bool isvol, int ordr)
{
- return getMem();
+ // Set flags
+ isVol = isvol;
+ order = ordr;
+ hasVals = true;
+ // Update values
+ size_t ngp = getCache(order).getVGauss().getN();
+ Js.resize(ngp);
+ iJs.resize(ngp);
+ detJs.resize(ngp);
+ this->update();
}
-
/**
- * @brief Return cache holding shape functions and Gauss points for all Tri3 elements
+ * @brief Initialize precomputed gradients
*/
-Cache &Tri3::getVCache() const
+void Tri3::initGradients()
{
- return getCache(order);
+ // Set flags
+ hasGrads = true;
+ // Update values
+ size_t ngp = getCache(order).getVGauss().getN();
+ gradJs.resize(ngp);
+ gradDetJs.resize(ngp);
+ this->update();
+}
+/**
+ * @brief Update precomputed values and gradients
+ */
+void Tri3::update()
+{
+ // Update values
+ if (hasVals)
+ {
+ // Jacobian and normal
+ // Tri3 is not of same dimension as Problem (surface element)
+ if (!isVol)
+ {
+ this->buildSurfaceJ();
+ this->computeNormal();
+ }
+ // Tri3 is of same dimension as Problem (volume element)
+ else
+ this->buildJ();
+ // Volume and CG
+ this->computeV();
+ this->computeCg();
+ }
+ // Update gradients
+ if (hasGrads)
+ {
+ // Jacobian, volume and normal
+ // Tri3 is not of same dimension as Problem (surface element)
+ if (!isVol)
+ {
+ this->buildGradientSurfaceJ();
+ this->buildGradientV(3);
+ this->computeGradientNormal();
+ }
+ // Tri3 is of same dimension as Problem (volume element)
+ else
+ {
+ this->buildGradientJ();
+ this->buildGradientV(2);
+ }
+ }
}
/**
- * @brief Compute Jacobian determinant at Gauss point k
+ * @brief Compute Jacobian determinant
*
* detJ = norm(x_i dN_i,xi X x_i dN_i,eta)
*/
-double Tri3::buildJ(size_t k) const
+void Tri3::buildSurfaceJ()
{
- std::vector<Eigen::Vector3d> ts = this->buildTangents(k);
- return ts[0].cross(ts[1]).norm();
+ for (size_t k = 0; k < getCache(order).getVGauss().getN(); ++k)
+ {
+ std::vector<Eigen::Vector3d> ts = this->buildTangents(k);
+ detJs[k] = ts[0].cross(ts[1]).norm();
+ }
}
/**
- * @brief Compute Jacobian matrix, inverse and determinant at Gauss point k
+ * @brief Compute Jacobian matrix, inverse and determinant
*
* J_ij(2,2) = diN_n * xj_n (n is the node)
*/
-double Tri3::buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const
+void Tri3::buildJ()
{
- Eigen::MatrixXd const &dff = getCache(order).getDsf(k);
+ CacheTri3 &cache = getCache(order);
+ for (size_t k = 0; k < cache.getVGauss().getN(); ++k)
+ {
+ Eigen::MatrixXd const &dff = getCache(order).getDsf(k);
- Eigen::Matrix2d JJ = Eigen::Matrix2d::Zero(); // temporary fixed-size matrix to efficiently compute the inverse
- size_t i = 0;
- for (auto it = nodes.begin(); it != nodes.end(); ++it, ++i)
- JJ += dff.col(i) * (*it)->pos.transpose().block<1, 2>(0, 0);
- J = JJ;
- iJ = JJ.inverse();
- return JJ.determinant();
+ Eigen::Matrix2d JJ = Eigen::Matrix2d::Zero(); // temporary fixed-size matrix to efficiently compute the inverse
+ size_t i = 0;
+ for (auto it = nodes.begin(); it != nodes.end(); ++it, ++i)
+ JJ += dff.col(i) * (*it)->pos.transpose().block<1, 2>(0, 0);
+ Js[k] = JJ;
+ iJs[k] = JJ.inverse();
+ detJs[k] = JJ.determinant();
+ this->checkJac(k);
+ }
}
/**
- * @brief Compute Jacobian determinant gradient with respect to each node coordinate, at Gauss point k
+ * @brief Compute Jacobian determinant gradient with respect to each node coordinate
*
* dDetJ[l = i * 3 + j] = detJ_l * (dt0 \cross t1 + t0 \cross dt1)_l;
*/
-void Tri3::buildGradientJ(size_t k, std::vector<double> &dDetJ) const
+void Tri3::buildGradientSurfaceJ()
{
- Eigen::MatrixXd const &dsf = getCache(order).getDsf(k);
-
- // Tangent vectors and determinant
- std::vector<Eigen::Vector3d> ts = this->buildTangents(k);
- Eigen::Vector3d detJ = ts[0].cross(ts[1]).normalized();
-
- // Jacobian gradient
- dDetJ.resize(3 * 3); // 3 nodes and 3 dimensions (x, y, z)
- for (size_t i = 0; i < 3; ++i)
+ CacheTri3 &cache = getCache(order);
+ for (size_t k = 0; k < getCache(order).getVGauss().getN(); ++k)
{
- for (size_t j = 0; j < 3; ++j)
+ // Shape functions, tangent vectors and determinant
+ Eigen::MatrixXd const &dsf = cache.getDsf(k);
+ std::vector<Eigen::Vector3d> ts = this->buildTangents(k);
+ Eigen::Vector3d detJ = ts[0].cross(ts[1]).normalized();
+
+ // Jacobian gradient
+ gradDetJs[k].resize(3 * 3); // 3 nodes and 3 dimensions (x, y, z)
+ for (size_t i = 0; i < 3; ++i)
{
- Eigen::Vector3d dt0 = Eigen::Vector3d::Zero(), dt1 = Eigen::Vector3d::Zero();
- dt0(j) = dsf(0, i);
- dt1(j) = dsf(1, i); // retain sf gradients associated to ith node for dimension j
- dDetJ[i * 3 + j] = detJ.dot(dt0.cross(ts[1]) + ts[0].cross(dt1));
+ for (size_t j = 0; j < 3; ++j)
+ {
+ Eigen::Vector3d dt0 = Eigen::Vector3d::Zero(), dt1 = Eigen::Vector3d::Zero();
+ dt0(j) = dsf(0, i);
+ dt1(j) = dsf(1, i); // retain sf gradients associated to ith node for dimension j
+ gradDetJs[k][i * 3 + j] = detJ.dot(dt0.cross(ts[1]) + ts[0].cross(dt1));
+ }
}
}
}
/**
- * @brief Compute Jacobian matrix gradient with respect to each node coordinate, at Gauss point k
+ * @brief Compute Jacobian matrix gradient with respect to each node coordinate
*
* dJ[l = i * 2 + j](2,2) = dN_i(2,1) * dxy_i(1,2)/dx_l (for node i and dimension j)
*/
-void Tri3::buildGradientJ(size_t k, std::vector<double> &dDetJ, std::vector<Eigen::MatrixXd> &dJ) const
+void Tri3::buildGradientJ()
{
- MemTri3 &mem = getMem();
- Eigen::MatrixXd const &dsf = getCache(order).getDsf(k);
-
- // Jacobian
- dJ.resize(3 * 2); // 3 nodes and 2 dimensions (x, y)
- for (size_t i = 0; i < 3; ++i)
+ CacheTri3 &cache = getCache(order);
+ for (size_t k = 0; k < cache.getVGauss().getN(); ++k)
{
- for (size_t j = 0; j < 2; ++j)
+ Eigen::MatrixXd const &dsf = cache.getDsf(k);
+
+ // Jacobian
+ gradJs[k].resize(3 * 2); // 3 nodes and 2 dimensions (x, y)
+ for (size_t i = 0; i < 3; ++i)
{
- dJ[i * 2 + j] = Eigen::Matrix2d::Zero();
- dJ[i * 2 + j].col(j) = dsf.col(i); // retain sf gradients associated to ith node for dimension j
+ for (size_t j = 0; j < 2; ++j)
+ {
+ gradJs[k][i * 2 + j] = Eigen::Matrix2d::Zero();
+ gradJs[k][i * 2 + j].col(j) = dsf.col(i); // retain sf gradients associated to ith node for dimension j
+ }
}
+ // Determinant
+ double detJ = this->getDetJ(k);
+ Eigen::MatrixXd const &iJ = this->getJinv(k);
+ gradDetJs[k].resize(3 * 2);
+ for (size_t i = 0; i < 3 * 2; ++i)
+ gradDetJs[k][i] = detJ * (iJ * gradJs[k][i]).trace();
}
- // Determinant
- double detJ = mem.getDetJ(k);
- Eigen::MatrixXd const &iJ = mem.getJinv(k);
- dDetJ.resize(3 * 2);
- for (size_t i = 0; i < dDetJ.size(); ++i)
- dDetJ[i] = detJ * (iJ * dJ[i]).trace();
}
/**
* @brief Return the element surface
*/
-double Tri3::computeV() const
+void Tri3::computeV()
{
CacheTri3 &cache = getCache(order);
GaussTri3 &gauss = cache.gauss;
- MemTri3 &mem = getMem();
// Gauss integration
- double V = 0.0;
+ vol = 0.0;
for (size_t k = 0; k < gauss.getN(); ++k)
- V += gauss.getW(k) * mem.getDetJ(k);
- return V;
+ vol += gauss.getW(k) * this->getDetJ(k);
}
/**
* @brief Compute volume gradient with respect to each node coordinate
*
* dV[i] = w_k * dDetJ[i]_k
*/
-std::vector<double> Tri3::buildGradientV(int dim) const
+void Tri3::buildGradientV(int dim)
{
GaussTri3 &gauss = getCache(order).gauss;
- MemTri3 &mem = getMem();
// Gradient of volume
- std::vector<double> dVol(3 * dim, 0.); // 3 nodes
+ gradVol.resize(3 * dim, 0.); // 3 nodes
for (size_t k = 0; k < gauss.getN(); ++k)
{
- double w = gauss.getW(k); // GP weight
- std::vector<double> dDetJ = mem.getGradDetJ(k); // gradient of Jacobian determinant
+ double w = gauss.getW(k); // GP weight
+ std::vector<double> dDetJ = this->getGradDetJ(k); // gradient of Jacobian determinant
for (size_t i = 0; i < dDetJ.size(); ++i)
- dVol[i] += w * dDetJ[i];
+ gradVol[i] += w * dDetJ[i];
}
- return dVol;
}
/**
@@ -201,16 +267,16 @@ std::vector<Eigen::Vector3d> Tri3::buildTangents(size_t k) const
/**
* @brief Compute element unit normal vector
*/
-Eigen::Vector3d Tri3::normal() const
+void Tri3::computeNormal()
{
std::vector<Eigen::Vector3d> ts = this->computeTangents();
- return ts[0].cross(ts[1]).normalized();
+ normal = ts[0].cross(ts[1]).normalized();
}
/**
* @brief Compute gradient of element unit normal vector with respect to each node coordinate
*/
-std::vector<Eigen::Vector3d> Tri3::gradientNormal() const
+void Tri3::computeGradientNormal()
{
// Compute normal and norm
std::vector<Eigen::Vector3d> ts = this->computeTangents();
@@ -219,17 +285,16 @@ std::vector<Eigen::Vector3d> Tri3::gradientNormal() const
// Compute gradients
Eigen::Matrix3d inn = (Eigen::Matrix3d::Identity() - n * n.transpose() / (nn * nn)) / nn; // gradient of n/|n|
- std::vector<Eigen::Vector3d> dn(3 * 3); // 3 nodes and 3 dimensions
+ gradNormal.resize(3 * 3); // 3 nodes and 3 dimensions
for (size_t i = 0; i < 3; ++i)
{
for (size_t j = 0; j < 3; ++j)
{
Eigen::Matrix3d dx = Eigen::Matrix3d::Zero();
- dx(j, i) = 1; // dpos_i / dx_j
- dn[i * 3 + j] = inn * ((dx.col(1) - dx.col(0)).cross(ts[1]) + (ts[0]).cross(dx.col(2) - dx.col(0))); // outer gradient times gradient of cross product
+ dx(j, i) = 1; // dpos_i / dx_j
+ gradNormal[i * 3 + j] = inn * ((dx.col(1) - dx.col(0)).cross(ts[1]) + (ts[0]).cross(dx.col(2) - dx.col(0))); // outer gradient times gradient of cross product
}
}
- return dn;
}
/**
@@ -242,5 +307,5 @@ Eigen::VectorXd Tri3::computeGradient(std::vector<double> const &u, size_t k) co
for (size_t i = 0; i < nodes.size(); ++i)
ue(i) = u[nodes[i]->row];
// Gradient in global axis: inv(J)_ij * d_jN_i * ue_i
- return getMem().getJinv(k) * getCache(order).getDsf(k) * ue;
+ return this->getJinv(k) * getCache(order).getDsf(k) * ue;
}
diff --git a/tbox/src/wTri3.h b/tbox/src/wTri3.h
index bb82eab1..c042c4c5 100644
--- a/tbox/src/wTri3.h
+++ b/tbox/src/wTri3.h
@@ -32,35 +32,31 @@ namespace tbox
class TBOX_API Tri3 : public Element
{
#ifndef SWIG
- friend class MemTri3;
-
private:
- MemTri3 *pmem; ///< private memory of precalculated values
- inline MemTri3 &getMem() const;
inline static CacheTri3 &getCache(int order);
protected:
- virtual double buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const override;
- virtual double buildJ(size_t k) const override;
- virtual void buildGradientJ(size_t k, std::vector<double> &dDetJ, std::vector<Eigen::MatrixXd> &dJ) const override;
- virtual void buildGradientJ(size_t k, std::vector<double> &dDetJ) const override;
-
-public: // @todo
- virtual double computeV() const override;
- virtual std::vector<double> buildGradientV(int dim) const override;
- virtual Eigen::Vector3d normal() const override;
- virtual std::vector<Eigen::Vector3d> gradientNormal() const override;
+ virtual void buildJ() override;
+ virtual void buildSurfaceJ() override;
+ virtual void buildGradientJ() override;
+ virtual void buildGradientSurfaceJ() override;
+ virtual void computeV() override;
+ virtual void buildGradientV(int dim) override;
+ virtual void computeNormal() override;
+ virtual void computeGradientNormal() override;
virtual std::vector<Eigen::Vector3d> computeTangents() const override;
virtual std::vector<Eigen::Vector3d> buildTangents(size_t k) const override;
#endif
public:
Tri3(int n, Tag *_ptag, Tag *_etag, std::vector<Tag *> &_parts, std::vector<Node *> &nods);
- virtual ~Tri3();
+ virtual ~Tri3() {}
virtual ELTYPE type() const override { return ELTYPE::TRI3; }
#ifndef SWIG
- virtual Mem &getVMem() const override;
+ virtual void initValues(bool isvol, int ordr = 2) override;
+ virtual void initGradients() override;
+ virtual void update() override;
virtual Cache &getVCache() const override;
virtual Eigen::VectorXd computeGradient(std::vector<double> const &u, size_t k) const override;
#endif
diff --git a/tbox/src/wTri3.inl b/tbox/src/wTri3.inl
index 23df4f1e..57f12736 100644
--- a/tbox/src/wTri3.inl
+++ b/tbox/src/wTri3.inl
@@ -14,14 +14,6 @@
* limitations under the License.
*/
-/**
- * @brief Return private memory for Tri3 private methods
- */
-inline MemTri3 &Tri3::getMem() const
-{
- return *pmem;
-}
-
/**
* @brief Return cache for Tri3 private methods
*/
diff --git a/tbox/src/wUnitTest.cpp b/tbox/src/wUnitTest.cpp
index e3f50dae..7b886c29 100644
--- a/tbox/src/wUnitTest.cpp
+++ b/tbox/src/wUnitTest.cpp
@@ -22,7 +22,6 @@
#include "wGmshExport.h"
#include "wResults.h"
#include "wMshDeform.h"
-#include "wMem.h"
#include <Eigen/Dense>
using namespace tbox;
@@ -84,7 +83,7 @@ void UnitTest::mshDeform(MshDeform *mshDef, int d)
// Initialize memory (volume element only)
for (auto e : mshDef->msh->elems)
if ((d == 2 && e->type() == ELTYPE::TRI3) || (d == 3 && e->type() == ELTYPE::TETRA4))
- e->getVMem().update(true);
+ e->initValues(true);
// Call deform
mshDef->deform();
}
diff --git a/waves/src/wProblem.cpp b/waves/src/wProblem.cpp
index b7a4a891..82907475 100644
--- a/waves/src/wProblem.cpp
+++ b/waves/src/wProblem.cpp
@@ -20,7 +20,7 @@
#include "wMedium.h"
#include "wBoundary.h"
#include "wElement.h"
-#include "wMem.h"
+
#include "wTag.h"
using namespace waves;
@@ -37,19 +37,18 @@ Problem::~Problem()
}
/**
- * @brief Update the elements memory (Jacobian)
- * @authors Adrien Crovato
+ * @brief Initialize the elements precomputed values
*/
-void Problem::updateMem()
+void Problem::initElems()
{
// Update volume Jacobian
for (auto vol : media)
for (auto e : vol->tag->elems)
- e->getVMem().update(true);
+ e->initValues(true);
// Update surface Jacobian (Neumann B.C.)
for (auto surf : bnds)
for (auto e : surf->tag->elems)
- e->getVMem().update(false);
+ e->initValues(false);
}
void Problem::write(std::ostream &out) const
diff --git a/waves/src/wProblem.h b/waves/src/wProblem.h
index aaf5b963..0f9e9b04 100644
--- a/waves/src/wProblem.h
+++ b/waves/src/wProblem.h
@@ -49,7 +49,7 @@ public:
void add(std::shared_ptr<Boundary> b) { bnds.push_back(b); }
#ifndef SWIG
- void updateMem();
+ void initElems();
virtual void write(std::ostream &out) const override;
#endif
};
diff --git a/waves/src/wTimeIntegration.cpp b/waves/src/wTimeIntegration.cpp
index d9586273..9d3c476c 100644
--- a/waves/src/wTimeIntegration.cpp
+++ b/waves/src/wTimeIntegration.cpp
@@ -52,7 +52,7 @@ TimeIntegration::TimeIntegration(std::shared_ptr<Problem> _pbl) : pbl(_pbl)
stopit = false;
// Update element memory
- pbl->updateMem();
+ pbl->initElems();
}
void TimeIntegration::start(MshExport *mshWriter)
@@ -101,7 +101,7 @@ void TimeIntegration::buildS(Eigen::SparseMatrix<double, Eigen::RowMajor> &S2)
for (auto bnd : pbl->bnds)
{
std::cout << "\tprocessing " << *bnd << '\n';
- std::cout << "normal=" << static_cast<Quad4 *>(bnd->tag->elems[0])->normal() << '\n';
+ std::cout << "normal=" << static_cast<Quad4 *>(bnd->tag->elems[0])->getNormal() << '\n';
tbb::parallel_do(bnd->tag->elems.begin(), bnd->tag->elems.end(),
[&](Element *e) {
if (e->type() != ELTYPE::QUAD4)
diff --git a/waves/src/wWaveTerm.cpp b/waves/src/wWaveTerm.cpp
index e39582ec..8cb036fc 100644
--- a/waves/src/wWaveTerm.cpp
+++ b/waves/src/wWaveTerm.cpp
@@ -19,7 +19,6 @@
#include "wElement.h"
#include "wCache.h"
#include "wGauss.h"
-#include "wMem.h"
using namespace waves;
using namespace tbox;
@@ -32,18 +31,17 @@ Eigen::MatrixXd WaveTerm::buildK(Element const &e, std::vector<double> const &u)
// Get precomputed values
Cache &cache = e.getVCache();
Gauss &gauss = cache.getVGauss();
- Mem &mem = e.getVMem();
// Stiffness matrix
Eigen::MatrixXd K = Eigen::MatrixXd::Zero(e.nodes.size(), e.nodes.size());
for (size_t k = 0; k < gauss.getN(); ++k)
{
// Jacobian inverse and shape functions
- Eigen::Matrix3d const &J = mem.getJinv(k);
+ Eigen::Matrix3d const &J = e.getJinv(k);
Eigen::MatrixXd const &dff = cache.getDsf(k);
// Elementary stiffness matrix
- K += (J * dff).transpose() * J * dff * gauss.getW(k) * mem.getDetJ(k);
+ K += (J * dff).transpose() * J * dff * gauss.getW(k) * e.getDetJ(k);
}
return K;
}
@@ -56,7 +54,6 @@ Eigen::MatrixXd WaveTerm::buildM(Element const &e)
// Get precomputed values
Cache &cache = e.getVCache();
Gauss &gauss = cache.getVGauss();
- Mem &mem = e.getVMem();
// Mass matrix
Eigen::MatrixXd M = Eigen::MatrixXd::Zero(e.nodes.size(), e.nodes.size());
@@ -64,7 +61,7 @@ Eigen::MatrixXd WaveTerm::buildM(Element const &e)
{
// Shape functions
Eigen::VectorXd const &ff = cache.getSf(k);
- M += ff * ff.transpose() * gauss.getW(k) * mem.getDetJ(k);
+ M += ff * ff.transpose() * gauss.getW(k) * e.getDetJ(k);
}
return M;
}
--
GitLab