From 8ad2c792a4315b40fdbfdde6d65b9a05a4d50639 Mon Sep 17 00:00:00 2001
From: acrovato <170-acrovato@users.noreply.gitlab.uliege.be>
Date: Sat, 21 Mar 2020 21:35:47 +0100
Subject: [PATCH] Rewritten Element methods. @K.Liegeois work remains.
---
mirrors/src/wSolver.cpp | 37 +++---
tbox/src/wElement.cpp | 17 +--
tbox/src/wElement.h | 7 +-
tbox/src/wHex8.cpp | 158 ++++++++++++------------
tbox/src/wHex8.h | 6 +-
tbox/src/wLine2.cpp | 55 ++++-----
tbox/src/wMshDeform.cpp | 2 +-
tbox/src/wQuad4.cpp | 48 +++++---
tbox/src/wTetra4.cpp | 216 +++++++++++----------------------
tbox/src/wTetra4.h | 7 +-
tbox/src/wTri3.cpp | 260 ++++++++++++++++++----------------------
tbox/src/wTri3.h | 2 +-
12 files changed, 355 insertions(+), 460 deletions(-)
diff --git a/mirrors/src/wSolver.cpp b/mirrors/src/wSolver.cpp
index 17c91fe4..d783dc87 100644
--- a/mirrors/src/wSolver.cpp
+++ b/mirrors/src/wSolver.cpp
@@ -161,23 +161,20 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
double E = m->E;
double nu = m->nu;
- double lambda = (nu * E) / ((1 + nu) * (1 - 2 * nu));
- double G = E / (2 * (1 + nu));
-
- Eigen::MatrixXd H = Eigen::MatrixXd::Zero(9, 9);
-
- for (auto i = 0; i < 3; ++i)
- for (auto j = 0; j < 3; ++j)
- for (auto k = 0; k < 3; ++k)
- for (auto l = 0; l < 3; ++l)
- {
- if (i == j && k == l)
- H(i * 3 + j, k * 3 + l) = lambda;
- if (i == k && j == l)
- H(i * 3 + j, k * 3 + l) = H(i * 3 + j, k * 3 + l) + G;
- if (i == l && j == k)
- H(i * 3 + j, k * 3 + l) = H(i * 3 + j, k * 3 + l) + G;
- }
+ Eigen::MatrixXd H = Eigen::MatrixXd::Zero(6, 6);
+ H(0, 0) = 1 - nu;
+ H(1, 1) = 1 - nu;
+ H(2, 2) = 1 - nu;
+ H(3, 3) = (1 - 2 * nu) / 2;
+ H(4, 4) = (1 - 2 * nu) / 2;
+ H(5, 5) = (1 - 2 * nu) / 2;
+ H(0, 1) = nu;
+ H(1, 0) = H(0, 1);
+ H(0, 2) = nu;
+ H(2, 0) = H(0, 2);
+ H(1, 2) = nu;
+ H(2, 1) = H(1, 2);
+ H *= E / ((1 + nu) * (1 - 2 * nu));
// K matrix
//=============================tbb==================================
@@ -192,7 +189,7 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
if (e->type() != ELTYPE::HEX8 && e->type() != ELTYPE::TETRA4)
return; //if the element is not a tetrahedron or a hexahedron, it is not treated
- e->buildK_mechanic(K_e, H);
+ e->buildK_mechanical(K_e, H);
for (size_t ii = 0; ii < e->nodes.size(); ++ii)
{
@@ -229,7 +226,7 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
if (e->type() != ELTYPE::HEX8 && e->type() != ELTYPE::TETRA4)
return; //if the element is not a tetrahedron or a hexahedron, it is not treated
- e->buildL_thermic(L_e, K_conductivity);
+ e->buildL_thermal(L_e, K_conductivity);
for (size_t ii = 0; ii < e->nodes.size(); ++ii)
{
@@ -271,7 +268,7 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
if (e->type() != ELTYPE::HEX8 && e->type() != ELTYPE::TETRA4)
return; //if the element is not a tetrahedron or a hexahedron, it is not treated
- e->buildS_thermomechanic(S_e, D);
+ e->buildS_thermomechanical(S_e, D);
for (size_t ii = 0; ii < e->nodes.size(); ++ii)
{
diff --git a/tbox/src/wElement.cpp b/tbox/src/wElement.cpp
index cabcb7a6..0105fcca 100644
--- a/tbox/src/wElement.cpp
+++ b/tbox/src/wElement.cpp
@@ -204,24 +204,19 @@ void Element::buildDK(Eigen::MatrixXd &K, std::vector<double> const &u, Fct0 con
throw std::runtime_error("Element::buildNK not implemented\n");
}
//---------------------
-void Element::buildK_linElastic(Eigen::MatrixXd &K, Eigen::MatrixXd const &H)
+void Element::buildK_mechanical(Eigen::MatrixXd &K_e, Eigen::MatrixXd const &H)
{
- throw std::runtime_error("Element::buildK_linElastic not implemented\n");
+ throw std::runtime_error("Element::buildK_mechanical not implemented\n");
}
-void Element::buildK_mechanic(Eigen::MatrixXd &K_e, Eigen::MatrixXd const &H)
+void Element::buildS_thermomechanical(Eigen::MatrixXd &S_e, Eigen::MatrixXd const &D)
{
- throw std::runtime_error("Element::buildK_mechanic not implemented\n");
+ throw std::runtime_error("Element::buildS_thermomechanical not implemented\n");
}
-void Element::buildS_thermomechanic(Eigen::MatrixXd &S_e, Eigen::MatrixXd const &D)
+void Element::buildL_thermal(Eigen::MatrixXd &L_e, Eigen::MatrixXd const &C)
{
- throw std::runtime_error("Element::buildS_thermomechanic not implemented\n");
-}
-
-void Element::buildL_thermic(Eigen::MatrixXd &L_e, Eigen::MatrixXd const &K_conductivity)
-{
- throw std::runtime_error("Element::buildL_thermic not implemented\n");
+ throw std::runtime_error("Element::buildL_thermial not implemented\n");
}
void Element::strain_stress(Eigen::MatrixXd &Strain, Eigen::MatrixXd &Stress, Eigen::MatrixXd const &H, std::vector<double> const &u)
diff --git a/tbox/src/wElement.h b/tbox/src/wElement.h
index c44114cb..51276571 100644
--- a/tbox/src/wElement.h
+++ b/tbox/src/wElement.h
@@ -92,10 +92,9 @@ public:
virtual void buildK(Eigen::MatrixXd &K, std::vector<double> const &u, Fct2 const &fct, bool fake);
virtual void buildDK(Eigen::MatrixXd &K, std::vector<double> const &u, Fct0 const &fct);
//---------------------
- virtual void buildK_linElastic(Eigen::MatrixXd &K, Eigen::MatrixXd const &H);
- virtual void buildK_mechanic(Eigen::MatrixXd &K_e, Eigen::MatrixXd const &H);
- virtual void buildS_thermomechanic(Eigen::MatrixXd &S_e, Eigen::MatrixXd const &D);
- virtual void buildL_thermic(Eigen::MatrixXd &L_e, Eigen::MatrixXd const &K_conductivity);
+ virtual void buildK_mechanical(Eigen::MatrixXd &K_e, Eigen::MatrixXd const &H);
+ virtual void buildS_thermomechanical(Eigen::MatrixXd &S_e, Eigen::MatrixXd const &D);
+ virtual void buildL_thermal(Eigen::MatrixXd &L_e, Eigen::MatrixXd const &C);
virtual void strain_stress(Eigen::MatrixXd &Strain, Eigen::MatrixXd &Stress, Eigen::MatrixXd const &H, std::vector<double> const &u);
virtual void iso_coordinates(Pt x, double eta, double dzeta);
virtual Pt normal();
diff --git a/tbox/src/wHex8.cpp b/tbox/src/wHex8.cpp
index d0597d61..2e33566d 100644
--- a/tbox/src/wHex8.cpp
+++ b/tbox/src/wHex8.cpp
@@ -36,7 +36,7 @@ Hex8::~Hex8()
delete pmem;
}
-/*
+/**
* @brief Private memory handling Jacobian for Hex8 Element
* @authors Romain Boman, Adrien Crovato
*/
@@ -49,7 +49,7 @@ MemHex8 &Hex8::getMem()
return *pmem;
}
-/*
+/**
* @brief Private cache handling shape functions and Gauss points for Hex8
* @authors Romain Boman, Adrien Crovato
*/
@@ -65,6 +65,8 @@ Cache &Hex8::getVCache()
/**
* @brief Compute Jacobian Matrix at Gauss point k
+ *
+ * J_ij(3,3) = diN_n * xj_n (n is the node)
* @authors Romain Boman
*/
double Hex8::buildJ(size_t k, Eigen::MatrixXd &J) const
@@ -81,83 +83,78 @@ double Hex8::buildJ(size_t k, Eigen::MatrixXd &J) const
return J.determinant();
}
+/**
+ * @brief Compute stiffness matrix for scalar field
+ *
+ * K_ij(8,8) = int f*dN_j*dN_i dV
+ * @authors Romain Boman
+ */
void Hex8::buildK(Eigen::MatrixXd &K, std::vector<double> const &u, Fct0 const &fct)
{
CacheHex8 &cache = getCache();
GaussHex8 &gauss = cache.gauss;
MemHex8 &mem = getMem();
- K = Eigen::MatrixXd::Zero(8, 8);
- Eigen::Vector3d factor1;
- Eigen::Vector3d factor2;
+ // Gauss integration
+ K = Eigen::Matrix<double, 8, 8>::Zero();
for (size_t k = 0; k < gauss.p.size(); ++k)
{
- // function and Jacobian
- double fk = fct.eval(*this, u, k);
- double detJ = mem.getDetJ(k);
+ // Jacobian inverse and shape functions
Eigen::Matrix3d const &J = mem.getJinv(k);
-
+ Eigen::Matrix<double, 3, 8> dN;
for (size_t i = 0; i < nodes.size(); ++i)
- {
- std::vector<double> &dffi = cache.dff[k][i];
- factor1 = J * Eigen::Vector3d(dffi.data());
+ dN.col(i) = Eigen::Vector3d(cache.dff[k][i].data());
- for (size_t j = 0; j < nodes.size(); ++j)
- {
- std::vector<double> &dffj = cache.dff[k][j];
- factor2 = J * Eigen::Vector3d(dffj.data());
- K(i, j) += factor1.dot(factor2) * gauss.w[k] * detJ * fk;
- }
- }
+ // Elementary stiffness matrix
+ K += (J * dN).transpose() * J * dN * fct.eval(*this, u, k) * gauss.w[k] * mem.getDetJ(k);
}
}
-void Hex8::buildK_mechanic(Eigen::MatrixXd &K_e, Eigen::MatrixXd const &H)
+/**
+ * @brief Compute linear elasticity stiffness matrix
+ *
+ * K_ij(24,24) = int eps_k H_ijkl eps_l dV = B(24,6)H(6,6)B(6,24)
+ * @authors Adrien Crovato
+ */
+void Hex8::buildK_mechanical(Eigen::MatrixXd &K, Eigen::MatrixXd const &H)
{
CacheHex8 &cache = getCache();
GaussHex8 &gauss = cache.gauss;
MemHex8 &mem = getMem();
- K_e = Eigen::MatrixXd::Zero(24, 24);
- Eigen::Vector3d factor1;
- Eigen::Vector3d factor2;
- Eigen::Vector3d factor3;
- Eigen::Matrix3d TMP1;
- for (size_t a = 0; a < gauss.p.size(); ++a)
+ K = Eigen::MatrixXd::Zero(24, 24);
+ for (size_t k = 0; k < gauss.p.size(); ++k)
{
- // calcul de la matrice jacobienne (au point de Gauss "ksi")
- double detJ = mem.getDetJ(a);
- Eigen::Matrix3d const &J = mem.getJinv(a);
-
- // calcul de K (au point de Gauss "ksi")
+ // Jacobian inverse
+ Eigen::Matrix3d const &invJ = mem.getJinv(k);
+ // Fill B matrix (shape functions)
+ Eigen::MatrixXd B = Eigen::MatrixXd::Zero(6, 3 * nodes.size());
for (size_t i = 0; i < nodes.size(); ++i)
{
- std::vector<double> &dffi = cache.dff[a][i];
- factor1 = J * Eigen::Vector3d(dffi.data());
-
- for (size_t j = i; j < nodes.size(); ++j)
- {
- std::vector<double> &dffj = cache.dff[a][j];
- factor2 = J * Eigen::Vector3d(dffj.data());
-
- for (size_t k = 0; k < 3; ++k)
- for (size_t m = 0; m < 3; ++m)
- {
- for (size_t l = 0; l < 3; ++l)
- for (size_t n = 0; n < 3; ++n)
- TMP1(l, n) = H(k * 3 + l, m * 3 + n);
- factor3 = TMP1 * factor2;
- K_e(i * 3 + k, j * 3 + m) += factor1.dot(factor3) * gauss.w[a] * detJ;
- }
- }
+ // Compute [Ni,x Ni,y, Ni_z]
+ std::vector<double> &dffi = cache.dff[k][i];
+ Eigen::Vector3d dN = invJ * Eigen::Vector3d(dffi.data());
+ B(0, 3 * i) = dN(0);
+ B(1, 3 * i + 1) = dN(1);
+ B(2, 3 * i + 2) = dN(2);
+ B(3, 3 * i + 1) = dN(2);
+ B(3, 3 * i + 2) = dN(1);
+ B(4, 3 * i) = dN(2);
+ B(4, 3 * i + 2) = dN(0);
+ B(5, 3 * i) = dN(1);
+ B(5, 3 * i + 1) = dN(0);
}
+
+ // Compute stiffness matrix
+ K += B.transpose() * H * B * gauss.w[k] * mem.getDetJ(k);
}
- for (size_t i = 0; i < 3 * nodes.size(); ++i)
- for (size_t j = i; j < 3 * nodes.size(); ++j)
- K_e(j, i) = K_e(i, j);
}
-void Hex8::buildS_thermomechanic(Eigen::MatrixXd &S_e, Eigen::MatrixXd const &D)
+/**
+ * @brief Compute thermomecanical matrix
+ * @authors Kim Liegeois
+ */
+void Hex8::buildS_thermomechanical(Eigen::MatrixXd &S_e, Eigen::MatrixXd const &D)
{
CacheHex8 &cache = getCache();
GaussHex8 &gauss = cache.gauss;
@@ -186,33 +183,33 @@ void Hex8::buildS_thermomechanic(Eigen::MatrixXd &S_e, Eigen::MatrixXd const &D)
}
}
-void Hex8::buildL_thermic(Eigen::MatrixXd &L_e, Eigen::MatrixXd const &K_conductivity)
+/**
+ * @brief Compute thermal matrix
+ * @authors Kim Liegeois
+ */
+void Hex8::buildL_thermal(Eigen::MatrixXd &L_e, Eigen::MatrixXd const &C)
{
CacheHex8 &cache = getCache();
GaussHex8 &gauss = cache.gauss;
MemHex8 &mem = getMem();
- L_e = Eigen::MatrixXd::Zero(8, 8);
- Eigen::Vector3d factor1;
- Eigen::Vector3d factor2;
- Eigen::Vector3d factor3;
+ // Gauss integration
+ L_e = Eigen::Matrix<double, 8, 8>::Zero();
for (size_t a = 0; a < gauss.p.size(); ++a)
{
- // calcul de la matrice jacobienne (au point de Gauss "ksi")
+ // Jacobian
double detJ = mem.getDetJ(a);
Eigen::Matrix3d const &J = mem.getJinv(a);
- // calcul de K (au point de Gauss "ksi")
for (size_t i = 0; i < nodes.size(); ++i)
{
std::vector<double> &dffi = cache.dff[a][i];
- factor1 = J * Eigen::Vector3d(dffi.data());
- factor2 = K_conductivity * factor1;
+ Eigen::Vector3d factor1 = C * J * Eigen::Vector3d(dffi.data());
for (size_t j = i; j < nodes.size(); ++j)
{
std::vector<double> &dffj = cache.dff[a][j];
- factor3 = J * Eigen::Vector3d(dffj.data());
- L_e(i, j) += factor2.dot(factor3) * gauss.w[a] * detJ;
+ Eigen::Vector3d factor2 = J * Eigen::Vector3d(dffj.data());
+ L_e(i, j) += factor1.dot(factor2) * gauss.w[a] * detJ;
}
}
}
@@ -264,22 +261,34 @@ void Hex8::strain_stress(Eigen::MatrixXd &Strain, Eigen::MatrixXd &Stress, Eigen
Stress(i, j) += H(i * 3 + j, k * 3 + l) * Strain(k, l);
}
+/**
+ * @brief Compute mass matrix
+ *
+ * M_ij(8,8) = int N_i*N_j dV
+ * @authors Romain Boman
+ */
void Hex8::buildM(Eigen::MatrixXd &M)
{
CacheHex8 &cache = getCache();
GaussHex8 &gauss = cache.gauss;
MemHex8 &mem = getMem();
- M = Eigen::MatrixXd::Zero(8, 8);
+ // Gauss integration
+ M = Eigen::Matrix<double, 8, 8>::Zero();
for (size_t k = 0; k < gauss.p.size(); ++k)
{
- // calcul de la matrice jacobienne (au point de Gauss "ksi")
+ // Jacobian and shape functions
double detJ = mem.getDetJ(k);
std::vector<double> &ff = cache.ff[k];
- M += Eigen::Matrix<double, 8, 1>(ff.data()) * Eigen::Matrix<double, 1, 8>(ff.data()) * gauss.w[k] * detJ; // M = M + {ff}^T [ff].wG.detJ
+
+ M += Eigen::Matrix<double, 8, 1>(ff.data()) * Eigen::Matrix<double, 1, 8>(ff.data()) * gauss.w[k] * detJ;
}
}
+/**
+ * @brief API for build methods
+ * @authors Romain Boman
+ */
void Hex8::build(MATTYPE type, Eigen::MatrixXd &A, std::vector<double> const &u, Fct0 const &fct)
{
switch (type)
@@ -298,8 +307,8 @@ void Hex8::build(MATTYPE type, Eigen::MatrixXd &A, std::vector<double> const &u,
}
/**
- * @brief Compute intergal of function
- * I = int fct dV
+ * @brief Compute intergal of scalar function
+ * I = int f dV
* @authors Adrien Crovato
*/
double Hex8::computeV(std::vector<double> const &u, Fct0 const &fct)
@@ -308,14 +317,10 @@ double Hex8::computeV(std::vector<double> const &u, Fct0 const &fct)
GaussHex8 &gauss = cache.gauss;
MemHex8 &mem = getMem();
+ // Gauss integration
double V = 0.0;
for (size_t k = 0; k < gauss.p.size(); ++k)
- {
- // function evaluation
- double fk = fct.eval(*this, u, k);
- double detJ = mem.getDetJ(k);
- V += fk * gauss.w[k] * detJ;
- }
+ V += fct.eval(*this, u, k) * gauss.w[k] * mem.getDetJ(k);
return V;
}
@@ -325,6 +330,5 @@ double Hex8::computeV(std::vector<double> const &u, Fct0 const &fct)
*/
double Hex8::computeV()
{
- MemHex8 &mem = getMem();
- return mem.getVol();
+ return getMem().getVol();
}
diff --git a/tbox/src/wHex8.h b/tbox/src/wHex8.h
index f3e86352..a3110597 100644
--- a/tbox/src/wHex8.h
+++ b/tbox/src/wHex8.h
@@ -39,9 +39,9 @@ class TBOX_API Hex8 : public Element
virtual void build(MATTYPE type, Eigen::MatrixXd &A, std::vector<double> const &u, Fct0 const &fct);
virtual void buildK(Eigen::MatrixXd &K, std::vector<double> const &u, Fct0 const &fct) override;
virtual void buildM(Eigen::MatrixXd &M) override;
- virtual void buildK_mechanic(Eigen::MatrixXd &K_e, Eigen::MatrixXd const &H) override;
- virtual void buildS_thermomechanic(Eigen::MatrixXd &S_e, Eigen::MatrixXd const &D) override;
- virtual void buildL_thermic(Eigen::MatrixXd &L_e, Eigen::MatrixXd const &K_conductivity) override;
+ virtual void buildK_mechanical(Eigen::MatrixXd &K_e, Eigen::MatrixXd const &H) override;
+ virtual void buildS_thermomechanical(Eigen::MatrixXd &S_e, Eigen::MatrixXd const &D) override;
+ virtual void buildL_thermal(Eigen::MatrixXd &L_e, Eigen::MatrixXd const &C) override;
virtual void strain_stress(Eigen::MatrixXd &Strain, Eigen::MatrixXd &Stress, Eigen::MatrixXd const &H, std::vector<double> const &u) override;
diff --git a/tbox/src/wLine2.cpp b/tbox/src/wLine2.cpp
index 1664ef6c..989ceb08 100644
--- a/tbox/src/wLine2.cpp
+++ b/tbox/src/wLine2.cpp
@@ -38,7 +38,7 @@ Line2::~Line2()
delete pmem;
}
-/*
+/**
* @brief Private memory handling Jacobian for Line2 Element
* @authors Romain Boman, Adrien Crovato
*/
@@ -51,7 +51,7 @@ MemLine2 &Line2::getMem()
return *pmem;
}
-/*
+/**
* @brief Private cache handling shape functions and Gauss points for Line2
* @authors Romain Boman, Adrien Crovato
*/
@@ -65,6 +65,10 @@ Cache &Line2::getVCache()
return getCache();
}
+/**
+ * @brief Compute element unit normal vector
+ * @authors Adrien Crovato
+ */
Pt Line2::normal()
{
Pt dx = (nodes[1]->pos - nodes[0]->pos).normalized();
@@ -73,6 +77,8 @@ Pt Line2::normal()
/**
* @brief Compute Jacombian determinant at Gauss point k
+ *
+ * detJ = norm(x_i * dN_i)
* @authors Romain Boman
*/
double Line2::buildJ(size_t k) const
@@ -82,10 +88,7 @@ double Line2::buildJ(size_t k) const
Eigen::Vector3d t0 = Eigen::Vector3d::Zero();
size_t i = 0;
for (auto it = nodes.begin(); it != nodes.end(); ++it, ++i)
- {
- std::vector<double> &dff = cache.dff[k][i];
- t0 += Eigen::Vector3d((*it)->pos.x.data()) * dff[0];
- }
+ t0 += Eigen::Vector3d((*it)->pos.x.data()) * cache.dff[k][i][0];
return t0.norm();
}
@@ -101,14 +104,10 @@ void Line2::builds(Eigen::VectorXd &s, std::vector<double> const &u, Fct0 const
GaussLine2 &gauss = cache.gauss;
MemLine2 &mem = getMem();
+ // Gauss integration
s = Eigen::Vector2d::Zero();
for (size_t k = 0; k < gauss.p.size(); ++k)
- {
- double sk = f.eval(*this, u, k);
- double detJ = mem.getDetJ(k);
-
- s += Eigen::Vector2d(cache.ff[k].data()) * sk * gauss.w[k] * detJ; // s = s + sk.wG.detJ
- }
+ s += Eigen::Vector2d(cache.ff[k].data()) * f.eval(*this, u, k) * gauss.w[k] * mem.getDetJ(k);
}
/**
@@ -123,31 +122,28 @@ void Line2::builds(Eigen::VectorXd &s, std::vector<double> const &u, Fct1 const
GaussLine2 &gauss = cache.gauss;
MemLine2 &mem = getMem();
+ // Gauss integration
s = Eigen::Vector2d::Zero();
for (size_t k = 0; k < gauss.p.size(); ++k)
- {
- double sk = f.eval(*this, u, k) * normal();
- double detJ = mem.getDetJ(k);
-
- s += Eigen::Vector2d(cache.ff[k].data()) * sk * gauss.w[k] * detJ; // s = s + sk.wG.detJ
- }
+ s += Eigen::Vector2d(cache.ff[k].data()) * (f.eval(*this, u, k) * normal()) * gauss.w[k] * mem.getDetJ(k);
}
/**
* @brief Compute square of gradient along element
+ *
+ * grad = (Delta_U / Delta_L)^2
* @authors Adrien Crovato
*/
double Line2::computeGrad2(std::vector<double> const &u, size_t k)
{
- double deltaU = u[nodes[1]->row] - u[nodes[0]->row];
- double deltaS = computeV();
- return deltaU / deltaS * deltaU / deltaS;
+ double dUdL = (u[nodes[1]->row] - u[nodes[0]->row]) / computeV(); // V = length of the element
+ return dUdL * dUdL;
}
/**
- * @brief Compute intergal of function
+ * @brief Compute intergal of scalar function
*
- * I = int fct dV
+ * I = int f dV
* @authors Adrien Crovato
*/
double Line2::computeV(std::vector<double> const &u, Fct0 const &fct)
@@ -156,23 +152,18 @@ double Line2::computeV(std::vector<double> const &u, Fct0 const &fct)
GaussLine2 &gauss = cache.gauss;
MemLine2 &mem = getMem();
+ // Gauss integration
double V = 0.0;
for (size_t k = 0; k < gauss.p.size(); ++k)
- {
- // evaluation de la fct a intégrer
- double fk = fct.eval(*this, u, k);
- double detJ = mem.getDetJ(k);
- V += fk * gauss.w[k] * detJ;
- }
+ V += fct.eval(*this, u, k) * gauss.w[k] * mem.getDetJ(k);
return V;
}
/**
- * @brief Return the element volume
+ * @brief Return the element length
* @authors Adrien Crovato
*/
double Line2::computeV()
{
- MemLine2 &mem = getMem();
- return mem.getVol();
+ return getMem().getVol();
}
diff --git a/tbox/src/wMshDeform.cpp b/tbox/src/wMshDeform.cpp
index 2b57c30a..e0b70a36 100644
--- a/tbox/src/wMshDeform.cpp
+++ b/tbox/src/wMshDeform.cpp
@@ -287,7 +287,7 @@ void MshDeform::deform()
// Elementary stiffness matrix
Eigen::MatrixXd Ke;
- e->buildK_linElastic(Ke, H);
+ e->buildK_mechanical(Ke, H);
// Assembly
// tbb::spin_mutex::scoped_lock lock(mutex);
diff --git a/tbox/src/wQuad4.cpp b/tbox/src/wQuad4.cpp
index 6129e1e3..6b350c5a 100644
--- a/tbox/src/wQuad4.cpp
+++ b/tbox/src/wQuad4.cpp
@@ -37,7 +37,7 @@ Quad4::~Quad4()
delete pmem;
}
-/*
+/**
* @brief Private memory handling Jacobian for Quad4 Element
* @authors Romain Boman, Adrien Crovato
*/
@@ -50,7 +50,7 @@ MemQuad4 &Quad4::getMem()
return *pmem;
}
-/*
+/**
* @brief Private cache handling shape functions and Gauss points for Quad4
* @authors Romain Boman, Adrien Crovato
*/
@@ -64,6 +64,10 @@ Cache &Quad4::getVCache()
return getCache();
}
+/**
+ * @brief Compute element unit normal vector
+ * @authors Romain Boman
+ */
Pt Quad4::normal()
{
Pt x1 = nodes[0]->pos;
@@ -76,6 +80,8 @@ Pt Quad4::normal()
/**
* @brief Compute Jacobian determinant at Gauss point k
+ *
+ * detJ = norm(x_i dN_i,xi X x_i dN_i,eta)
* @authors Romain Boman
*/
double Quad4::buildJ(size_t k) const
@@ -94,7 +100,9 @@ double Quad4::buildJ(size_t k) const
return t0.cross(t1).norm();
}
/**
- * @brief Compute Jacobian matrix at Gauss point k
+ * @brief Compute Jacobian matrix and determinant at Gauss point k
+ *
+ * J_ij(2,2) = diN_n * xj_n (n is the node)
* @authors Romain Boman
*/
double Quad4::buildJ(size_t k, Eigen::MatrixXd &J) const
@@ -111,23 +119,34 @@ double Quad4::buildJ(size_t k, Eigen::MatrixXd &J) const
return J.determinant();
}
+/**
+ * @brief Compute mass matrix
+ *
+ * S_ij(4,4) = int N_i*N_j dV (from Quad4)
+ * @authors Romain Boman
+ */
void Quad4::buildS(Eigen::MatrixXd &S)
{
CacheQuad4 &cache = getCache();
GaussQuad4 &gauss = cache.gauss;
MemQuad4 &mem = getMem();
+ // Gauss integration
S = Eigen::Matrix4d::Zero();
for (size_t k = 0; k < gauss.p.size(); ++k)
{
- // Jacobian
+ // Jacobian and shape function
double detJ = mem.getDetJ(k);
-
std::vector<double> &ff = cache.ff[k];
- S += Eigen::Vector4d(ff.data()) * Eigen::RowVector4d(ff.data()) * gauss.w[k] * detJ; // S = S + {ff}^T [ff].wG.detJ
+
+ S += Eigen::Vector4d(ff.data()) * Eigen::RowVector4d(ff.data()) * gauss.w[k] * detJ;
}
}
+/**
+ * @brief API for build methods
+ * @authors Romain Boman
+ */
void Quad4::build(MATTYPE type, Eigen::MatrixXd &A)
{
switch (type)
@@ -143,9 +162,9 @@ void Quad4::build(MATTYPE type, Eigen::MatrixXd &A)
}
/**
- * @brief Compute intergal of function
+ * @brief Compute intergal of scalar function
*
- * I = int fct dV
+ * I = int f dV
* @authors Romain Boman
*/
double Quad4::computeV(std::vector<double> const &u, Fct0 const &fct)
@@ -154,23 +173,18 @@ double Quad4::computeV(std::vector<double> const &u, Fct0 const &fct)
GaussQuad4 &gauss = cache.gauss;
MemQuad4 &mem = getMem();
+ // Gauss integration
double V = 0.0;
for (size_t k = 0; k < gauss.p.size(); ++k)
- {
- // Function evaluation
- double fk = fct.eval(*this, u, k);
- double detJ = mem.getDetJ(k);
- V += fk * gauss.w[k] * detJ;
- }
+ V += fct.eval(*this, u, k) * gauss.w[k] * mem.getDetJ(k);
return V;
}
/**
- * @brief Return the element volume
+ * @brief Return the element surface
* @authors Adrien Crovato
*/
double Quad4::computeV()
{
- MemQuad4 &mem = getMem();
- return mem.getVol();
+ return getMem().getVol();
}
\ No newline at end of file
diff --git a/tbox/src/wTetra4.cpp b/tbox/src/wTetra4.cpp
index c030915e..e2d996c1 100644
--- a/tbox/src/wTetra4.cpp
+++ b/tbox/src/wTetra4.cpp
@@ -36,7 +36,7 @@ Tetra4::~Tetra4()
delete pmem;
}
-/*
+/**
* @brief Private memory handling Jacobian for Tetra4 Element
* @authors Romain Boman, Adrien Crovato
*/
@@ -49,7 +49,7 @@ MemTetra4 &Tetra4::getMem()
return *pmem;
}
-/*
+/**
* @brief Private cache handling shape functions and Gauss points for Tetra4
* @authors Romain Boman, Adrien Crovato
*/
@@ -65,6 +65,8 @@ Cache &Tetra4::getVCache()
/**
* @brief Compute Jacobian Matrix at Gauss point k
+ *
+ * J_ij(3,3) = diN_n * xj_n (n is the node)
* @authors Romain Boman
*/
double Tetra4::buildJ(size_t k, Eigen::MatrixXd &J) const
@@ -83,6 +85,8 @@ double Tetra4::buildJ(size_t k, Eigen::MatrixXd &J) const
/**
* @brief Compute mass matrix
+ *
+ * M_ij(4,4) = int N_i*N_j dV
* @authors Romain Boman
*/
void Tetra4::buildM(Eigen::MatrixXd &M)
@@ -91,21 +95,22 @@ void Tetra4::buildM(Eigen::MatrixXd &M)
GaussTetra4 &gauss = cache.gauss;
MemTetra4 &mem = getMem();
+ // Gauss integration
M = Eigen::Matrix4d::Zero();
for (size_t k = 0; k < gauss.p.size(); ++k)
{
- // Jacobian
+ // Jacobian and shape functions
double detJ = mem.getDetJ(k);
-
std::vector<double> &ff = cache.ff[k];
- M += Eigen::Vector4d(ff.data()) * Eigen::RowVector4d(ff.data()) * gauss.w[k] * detJ; // M = M + {ff}^T [ff].wG.detJ
+
+ M += Eigen::Vector4d(ff.data()) * Eigen::RowVector4d(ff.data()) * gauss.w[k] * detJ;
}
}
/**
* @brief Compute stiffness matrix for scalar field
*
- * K_ij(4,4) = int fct*dN_j*dN_i dV
+ * K_ij(4,4) = int f*dN_j*dN_i dV
* @authors Romain Boman, Adrien Crovato
*/
void Tetra4::buildK(Eigen::MatrixXd &K, std::vector<double> const &u, Fct0 const &fct)
@@ -114,29 +119,18 @@ void Tetra4::buildK(Eigen::MatrixXd &K, std::vector<double> const &u, Fct0 const
GaussTetra4 &gauss = cache.gauss;
MemTetra4 &mem = getMem();
+ // Gauss integration
K = Eigen::Matrix4d::Zero();
- Eigen::Vector3d factor1;
- Eigen::Vector3d factor2;
for (size_t k = 0; k < gauss.p.size(); ++k)
{
- // Function and Jacobian evaluation
- double fk = fct.eval(*this, u, k);
- double detJ = mem.getDetJ(k);
+ // Jacobian inverse and shape functions
Eigen::Matrix3d const &J = mem.getJinv(k);
-
- // compute stiffness matrix for each GP
+ Eigen::Matrix<double, 3, 4> dN;
for (size_t i = 0; i < nodes.size(); ++i)
- {
- std::vector<double> &dffi = cache.dff[k][i];
- factor1 = J * Eigen::Vector3d(dffi.data());
+ dN.col(i) = Eigen::Vector3d(cache.dff[k][i].data());
- for (size_t j = 0; j < nodes.size(); ++j)
- {
- std::vector<double> &dffj = cache.dff[k][j];
- factor2 = J * Eigen::Vector3d(dffj.data());
- K(i, j) += factor1.dot(factor2) * gauss.w[k] * detJ * fk;
- }
- }
+ // Elementary stiffness matrix
+ K += (J * dN).transpose() * J * dN * fct.eval(*this, u, k) * gauss.w[k] * mem.getDetJ(k);
}
}
@@ -152,30 +146,20 @@ void Tetra4::buildDK(Eigen::MatrixXd &K, std::vector<double> const &u, Fct0 cons
GaussTetra4 &gauss = cache.gauss;
MemTetra4 &mem = getMem();
+ // Gauss integration
K = Eigen::Matrix4d::Zero();
- Eigen::Vector3d factor1;
- Eigen::Vector3d factor2;
for (size_t k = 0; k < gauss.p.size(); ++k)
{
- // "derivative" of function, Jacobian and gradient evaluation
- double dfk = fct.evalD(*this, u, k);
- double detJ = mem.getDetJ(k);
+ // Jacobian inverse, shape functions and gradient
Eigen::Matrix3d const &J = mem.getJinv(k);
+ Eigen::Matrix<double, 3, 4> dN;
+ for (size_t i = 0; i < nodes.size(); ++i)
+ dN.col(i) = Eigen::Vector3d(cache.dff[k][i].data());
Eigen::VectorXd gradu;
- computeGrad(u, gradu, k);
+ this->computeGrad(u, gradu, k);
- // computation of K
- for (size_t i = 0; i < nodes.size(); ++i)
- {
- std::vector<double> &dffi = cache.dff[k][i];
- factor1 = J * Eigen::Vector3d(dffi.data());
- for (size_t j = 0; j < nodes.size(); ++j)
- {
- std::vector<double> &dffj = cache.dff[k][j];
- factor2 = J * Eigen::Vector3d(dffj.data());
- K(i, j) += gradu.dot(factor1) * gradu.dot(factor2) * gauss.w[k] * detJ * dfk;
- }
- }
+ // Elementary stiffness matrix
+ K += (gradu.transpose() * J * dN).transpose() * gradu.transpose() * J * dN * fct.evalD(*this, u, k) * gauss.w[k] * mem.getDetJ(k);
}
}
@@ -191,181 +175,121 @@ void Tetra4::buildDs(Eigen::VectorXd &s, std::vector<double> const &u, Fct0 cons
GaussTetra4 &gauss = cache.gauss;
MemTetra4 &mem = getMem();
+ // Gauss integration
s = Eigen::Vector4d::Zero();
- Eigen::Vector3d factor1;
for (size_t k = 0; k < gauss.p.size(); ++k)
{
- // Gradient, function and Jacobian
+ // Shape functions and gradient
+ Eigen::Matrix<double, 3, 4> dN;
+ for (size_t i = 0; i < nodes.size(); ++i)
+ dN.col(i) = Eigen::Vector3d(cache.dff[k][i].data());
Eigen::VectorXd gradu;
this->computeGrad(u, gradu, k);
- double fk = f.eval(*this, u, k);
- double detJ = mem.getDetJ(k);
- Eigen::Matrix3d const &J = mem.getJinv(k);
- // s(i) = f * gradu * gradN(i)
- for (size_t i = 0; i < nodes.size(); ++i)
- {
- std::vector<double> &dffi = cache.dff[k][i];
- factor1 = J * Eigen::Vector3d(dffi.data());
- s(i) += gradu.dot(factor1) * gauss.w[k] * detJ * fk;
- }
+ // Elementary flux vector
+ s += gradu.transpose() * mem.getJinv(k) * dN * f.eval(*this, u, k) * gauss.w[k] * mem.getDetJ(k);
}
}
/**
* @brief Compute linear elasticity stiffness matrix
+ *
+ * K_ij(12,12) = int eps_k H_ijkl eps_l dV = B(12,6)H(6,6)B(6,12)
* @authors Adrien Crovato
- * @todo this is a clean duplicate of buildK_mechnical, should be merged
*/
-void Tetra4::buildK_linElastic(Eigen::MatrixXd &K, Eigen::MatrixXd const &H)
+void Tetra4::buildK_mechanical(Eigen::MatrixXd &K, Eigen::MatrixXd const &H)
{
CacheTetra4 &cache = getCache();
GaussTetra4 &gauss = cache.gauss;
MemTetra4 &mem = getMem();
K = Eigen::MatrixXd::Zero(12, 12);
- Eigen::Vector3d factori;
for (size_t k = 0; k < gauss.p.size(); ++k)
{
- // Jacobian
- double detJ = mem.getDetJ(k);
+ // Jacobian inverse
Eigen::Matrix3d const &invJ = mem.getJinv(k);
-
- // Fill B matrix
+ // Fill B matrix (shape functions)
Eigen::MatrixXd B = Eigen::MatrixXd::Zero(6, 3 * nodes.size());
for (size_t i = 0; i < nodes.size(); ++i)
{
// Compute [Ni,x Ni,y, Ni_z]
std::vector<double> &dffi = cache.dff[k][i];
- factori = invJ * Eigen::Vector3d(dffi.data());
- B(0, 3 * i) = factori(0);
- B(1, 3 * i + 1) = factori(1);
- B(2, 3 * i + 2) = factori(2);
- B(3, 3 * i + 1) = factori(2);
- B(3, 3 * i + 2) = factori(1);
- B(4, 3 * i) = factori(2);
- B(4, 3 * i + 2) = factori(0);
- B(5, 3 * i) = factori(1);
- B(5, 3 * i + 1) = factori(0);
+ Eigen::Vector3d dN = invJ * Eigen::Vector3d(dffi.data());
+ B(0, 3 * i) = dN(0);
+ B(1, 3 * i + 1) = dN(1);
+ B(2, 3 * i + 2) = dN(2);
+ B(3, 3 * i + 1) = dN(2);
+ B(3, 3 * i + 2) = dN(1);
+ B(4, 3 * i) = dN(2);
+ B(4, 3 * i + 2) = dN(0);
+ B(5, 3 * i) = dN(1);
+ B(5, 3 * i + 1) = dN(0);
}
// Compute stiffness matrix
- K += B.transpose() * H * B * detJ * gauss.w[k];
+ K += B.transpose() * H * B * gauss.w[k] * mem.getDetJ(k);
}
}
/* -------------------------------------------------------------------------- */
-/**
- * @brief Compute mechanical stiffness matrix
- * @authors Kim Liegeois
- */
-void Tetra4::buildK_mechanic(Eigen::MatrixXd &K_e, Eigen::MatrixXd const &H)
-{
- CacheTetra4 &cache = getCache();
- GaussTetra4 &gauss = cache.gauss;
- MemTetra4 &mem = getMem();
-
- K_e = Eigen::MatrixXd::Zero(12, 12);
- Eigen::Vector3d factor1;
- Eigen::Vector3d factor2;
- Eigen::Vector3d factor3;
- Eigen::Matrix3d TMP1;
- for (size_t a = 0; a < gauss.p.size(); ++a)
- {
- // calcul de la matrice jacobienne (au point de Gauss "ksi")
- double detJ = mem.getDetJ(a);
- Eigen::Matrix3d const &J = mem.getJinv(a);
-
- // calcul de K (au point de Gauss "ksi")
- for (size_t i = 0; i < nodes.size(); ++i)
- {
- std::vector<double> &dffi = cache.dff[a][i];
- factor1 = J * Eigen::Vector3d(dffi.data());
- for (size_t j = i; j < nodes.size(); ++j)
- {
- std::vector<double> &dffj = cache.dff[a][j];
- factor2 = J * Eigen::Vector3d(dffj.data());
- for (size_t k = 0; k < 3; ++k)
- for (size_t m = 0; m < 3; ++m)
- {
- for (size_t l = 0; l < 3; ++l)
- for (size_t n = 0; n < 3; ++n)
- TMP1(l, n) = H(k * 3 + l, m * 3 + n);
- factor3 = TMP1 * factor2;
- K_e(i * 3 + k, j * 3 + m) += factor1.dot(factor3) * gauss.w[a] * detJ;
- }
- }
- }
- }
- for (size_t i = 0; i < 3 * nodes.size(); ++i)
- for (size_t j = i; j < 3 * nodes.size(); ++j)
- K_e(j, i) = K_e(i, j);
-}
/**
* @brief Compute thermomechanical mass matrix
* @authors Kim Liegeois
*/
-void Tetra4::buildS_thermomechanic(Eigen::MatrixXd &S_e, Eigen::MatrixXd const &D)
+void Tetra4::buildS_thermomechanical(Eigen::MatrixXd &S_e, Eigen::MatrixXd const &D)
{
CacheTetra4 &cache = getCache();
GaussTetra4 &gauss = cache.gauss;
MemTetra4 &mem = getMem();
+ // Gauss integration
S_e = Eigen::MatrixXd::Zero(12, 4);
- Eigen::Vector3d factor1;
- Eigen::Vector3d factor2;
for (size_t a = 0; a < gauss.p.size(); ++a)
{
- // calcul de la matrice jacobienne (au point de Gauss "ksi")
+ // Jacobian and shape functions
double detJ = mem.getDetJ(a);
Eigen::Matrix3d const &J = mem.getJinv(a);
std::vector<double> &ff = cache.ff[a];
- // calcul de K (au point de Gauss "ksi")
for (size_t i = 0; i < nodes.size(); ++i)
{
std::vector<double> &dffi = cache.dff[a][i];
- factor1 = J * Eigen::Vector3d(dffi.data());
- factor2 = D * factor1;
+ Eigen::Vector3d factor = D * J * Eigen::Vector3d(dffi.data());
for (size_t j = 0; j < nodes.size(); ++j)
for (size_t k = 0; k < 3; ++k)
- S_e(i * 3 + k, j) += factor2[k] * ff[j] * gauss.w[a] * detJ;
+ S_e(i * 3 + k, j) += factor(k) * ff[j] * gauss.w[a] * detJ;
}
}
}
/**
- * @brief Compute thermic matrix
+ * @brief Compute thermal matrix
* @authors Kim Liegeois
*/
-void Tetra4::buildL_thermic(Eigen::MatrixXd &L_e, Eigen::MatrixXd const &K_conductivity)
+void Tetra4::buildL_thermal(Eigen::MatrixXd &L_e, Eigen::MatrixXd const &C)
{
CacheTetra4 &cache = getCache();
GaussTetra4 &gauss = cache.gauss;
MemTetra4 &mem = getMem();
+ // Gauss integration
L_e = Eigen::Matrix4d::Zero();
- Eigen::Vector3d factor1;
- Eigen::Vector3d factor2;
- Eigen::Vector3d factor3;
for (size_t a = 0; a < gauss.p.size(); ++a)
{
- // calcul de la matrice jacobienne (au point de Gauss "ksi")
+ // Jacobian
double detJ = mem.getDetJ(a);
Eigen::Matrix3d const &J = mem.getJinv(a);
- // calcul de K (au point de Gauss "ksi")
for (size_t i = 0; i < nodes.size(); ++i)
{
std::vector<double> &dffi = cache.dff[a][i];
- factor1 = J * Eigen::Vector3d(dffi.data());
- factor2 = K_conductivity * factor1;
+ Eigen::Vector3d factor1 = C * J * Eigen::Vector3d(dffi.data());
for (size_t j = i; j < nodes.size(); ++j)
{
std::vector<double> &dffj = cache.dff[a][j];
- factor3 = J * Eigen::Vector3d(dffj.data());
- L_e(i, j) += factor2.dot(factor3) * gauss.w[a] * detJ;
+ Eigen::Vector3d factor2 = J * Eigen::Vector3d(dffj.data());
+ L_e(i, j) += factor1.dot(factor2) * gauss.w[a] * detJ;
}
}
}
@@ -434,7 +358,7 @@ void Tetra4::computeGrad(std::vector<double> const &u, Eigen::VectorXd &grad, si
CacheTetra4 &cache = getCache();
MemTetra4 &mem = getMem();
- // local matrix of sf derivatives, solution vector and Jacobian
+ // local matrix of sf derivatives and solution vector
Eigen::MatrixXd gradN(3, nodes.size());
Eigen::Vector4d ue;
for (size_t i = 0; i < nodes.size(); ++i)
@@ -442,9 +366,8 @@ void Tetra4::computeGrad(std::vector<double> const &u, Eigen::VectorXd &grad, si
gradN.col(i) = Eigen::Vector3d(cache.dff[k][i].data());
ue(i) = u[nodes[i]->row];
}
- Eigen::Matrix3d const &Jinv = mem.getJinv(k);
// gradient in global axis: inv(J)_ij * d_jN_i * ue_i
- grad = Jinv * gradN * ue;
+ grad = mem.getJinv(k) * gradN * ue;
}
/**
@@ -459,9 +382,9 @@ double Tetra4::computeGrad2(std::vector<double> const &u, size_t k)
}
/**
- * @brief Compute intergal of function
+ * @brief Compute intergal of scalar function
*
- * I = int fct dV
+ * I = int f dV
* @authors Adrien Crovato
*/
double Tetra4::computeV(std::vector<double> const &u, Fct0 const &fct)
@@ -470,14 +393,10 @@ double Tetra4::computeV(std::vector<double> const &u, Fct0 const &fct)
GaussTetra4 &gauss = cache.gauss;
MemTetra4 &mem = getMem();
+ // Gauss integration
double V = 0.0;
for (size_t k = 0; k < gauss.p.size(); ++k)
- {
- // function evaluation
- double fk = fct.eval(*this, u, k);
- double detJ = mem.getDetJ(k);
- V += fk * gauss.w[k] * detJ;
- }
+ V += fct.eval(*this, u, k) * gauss.w[k] * mem.getDetJ(k);
return V;
}
@@ -487,6 +406,5 @@ double Tetra4::computeV(std::vector<double> const &u, Fct0 const &fct)
*/
double Tetra4::computeV()
{
- MemTetra4 &mem = getMem();
- return mem.getVol();
+ return getMem().getVol();
}
diff --git a/tbox/src/wTetra4.h b/tbox/src/wTetra4.h
index 17a10f90..3795dad6 100644
--- a/tbox/src/wTetra4.h
+++ b/tbox/src/wTetra4.h
@@ -40,12 +40,11 @@ class TBOX_API Tetra4 : public Element
virtual void buildK(Eigen::MatrixXd &K, std::vector<double> const &u, Fct0 const &fct) override;
virtual void buildDK(Eigen::MatrixXd &K, std::vector<double> const &u, Fct0 const &fct) override;
virtual void buildDs(Eigen::VectorXd &s, std::vector<double> const &u, Fct0 const &f) override;
- virtual void buildK_linElastic(Eigen::MatrixXd &K, Eigen::MatrixXd const &H) override;
+ virtual void buildK_mechanical(Eigen::MatrixXd &K, Eigen::MatrixXd const &H) override;
/* @todo Kim's work, will have to be cleaned at some point */
- virtual void buildK_mechanic(Eigen::MatrixXd &K_e, Eigen::MatrixXd const &H) override;
- virtual void buildS_thermomechanic(Eigen::MatrixXd &S_e, Eigen::MatrixXd const &D) override;
- virtual void buildL_thermic(Eigen::MatrixXd &L_e, Eigen::MatrixXd const &K_conductivity) override;
+ virtual void buildS_thermomechanical(Eigen::MatrixXd &S_e, Eigen::MatrixXd const &D) override;
+ virtual void buildL_thermal(Eigen::MatrixXd &L_e, Eigen::MatrixXd const &C) override;
virtual void strain_stress(Eigen::MatrixXd &Strain, Eigen::MatrixXd &Stress, Eigen::MatrixXd const &H, std::vector<double> const &u) override;
/* --- */
diff --git a/tbox/src/wTri3.cpp b/tbox/src/wTri3.cpp
index 961e2aa3..ee4391c2 100644
--- a/tbox/src/wTri3.cpp
+++ b/tbox/src/wTri3.cpp
@@ -39,7 +39,7 @@ Tri3::~Tri3()
delete pmem;
}
-/*
+/**
* @brief Private memory handling Jacobian for Tri3 Element
* @authors Romain Boman, Adrien Crovato
*/
@@ -52,7 +52,7 @@ MemTri3 &Tri3::getMem()
return *pmem;
}
-/*
+/**
* @brief Private cache handling shape functions and Gauss points for Tri3
* @authors Romain Boman, Adrien Crovato
*/
@@ -66,6 +66,10 @@ Cache &Tri3::getVCache()
return getCache();
}
+/**
+ * @brief Compute element unit normal vector
+ * @authors Adrien Crovato
+ */
Pt Tri3::normal()
{
Pt x1 = nodes[0]->pos;
@@ -77,7 +81,9 @@ Pt Tri3::normal()
/**
* @brief Compute Jacobian determinant at Gauss point k
- * @authors Romain Boman
+ *
+ * detJ = norm(x_i dN_i,xi X x_i dN_i,eta)
+ * @authors Romain Boman, Adrien Crovato
*/
double Tri3::buildJ(size_t k) const
{
@@ -95,8 +101,10 @@ double Tri3::buildJ(size_t k) const
return t0.cross(t1).norm();
}
/**
- * @brief Compute Jacobian matrix at Gauss point k
- * @authors Romain Boman
+ * @brief Compute Jacobian matrix and determinant at Gauss point k
+ *
+ * J_ij(2,2) = diN_n * xj_n (n is the node)
+ * @authors Romain Boman, Adrien Crovato
*/
double Tri3::buildJ(size_t k, Eigen::MatrixXd &J) const
{
@@ -115,8 +123,8 @@ double Tri3::buildJ(size_t k, Eigen::MatrixXd &J) const
/**
* @brief Compute stiffness matrix for scalar field
*
- * K_ij(3,3) = int fct*dN_j*dN_i dV
- * @authors Romain Boman
+ * K_ij(3,3) = int f*dN_j*dN_i dV
+ * @authors Romain Boman, Adrien Crovato
*/
void Tri3::buildK(Eigen::MatrixXd &K, std::vector<double> const &u, Fct0 const &fct)
{
@@ -124,38 +132,27 @@ void Tri3::buildK(Eigen::MatrixXd &K, std::vector<double> const &u, Fct0 const &
GaussTri3 &gauss = cache.gauss;
MemTri3 &mem = getMem();
+ // Gauss integration
K = Eigen::Matrix3d::Zero();
- Eigen::Vector2d factor1;
- Eigen::Vector2d factor2;
for (size_t k = 0; k < gauss.p.size(); ++k)
{
- // Function and Jacobian inverse
- double fk = fct.eval(*this, u, k);
- double detJ = mem.getDetJ(k);
+ // Jacobian inverse and shape functions
Eigen::Matrix2d const &J = mem.getJinv(k);
-
- // Elementary stiffness matrix
+ Eigen::Matrix<double, 2, 3> dN;
for (size_t i = 0; i < nodes.size(); ++i)
- {
- std::vector<double> &dffi = cache.dff[k][i];
- factor1 = J * Eigen::Vector2d(dffi.data());
-
- for (size_t j = 0; j < nodes.size(); ++j)
- {
- std::vector<double> &dffj = cache.dff[k][j];
- factor2 = J * Eigen::Vector2d(dffj.data());
+ dN.col(i) = Eigen::Vector2d(cache.dff[k][i].data());
- K(i, j) += factor1.dot(factor2) * gauss.w[k] * detJ * fk;
- }
- }
- // other way, for the future
- //Eigen::Matrix<double, 2, 3> dN;
- //for (size_t i = 0; i < nodes.size(); ++i)
- // dN.col(i) = Eigen::Vector2d(cache.dff[k][i].data());
- //K = (J * dN).transpose() * J * dN * gauss.w[k] * detJ * fk;
+ // Elementary stiffness matrix
+ K += (J * dN).transpose() * J * dN * fct.eval(*this, u, k) * gauss.w[k] * mem.getDetJ(k);
}
}
+/**
+ * @brief Compute stiffness matrix for scalar field
+ *
+ * K_ij(3,3) = int [f]*dN_j*dN_i dV
+ * @authors Romain Boman
+ */
void Tri3::buildK(Eigen::MatrixXd &K, std::vector<double> const &u, Fct2 const &fct, bool fake)
{
CacheTri3 &cache = getCache();
@@ -171,33 +168,21 @@ void Tri3::buildK(Eigen::MatrixXd &K, std::vector<double> const &u, Fct2 const &
}
else // true run - use MPI cached results
{
+ // Gauss integration
MemTri3 &mem = getMem();
K = Eigen::Matrix3d::Zero();
- Eigen::Vector2d factor0;
- Eigen::Vector2d factor1;
- Eigen::Vector2d factor2;
for (size_t k = 0; k < gauss.p.size(); ++k)
{
- // Function and Jacobian inverse
+ // Jacobian inverse, shape functions and function evaluation
+ Eigen::Matrix2d const &J = mem.getJinv(k);
+ Eigen::Matrix<double, 2, 3> dN;
+ for (size_t i = 0; i < nodes.size(); ++i)
+ dN.col(i) = Eigen::Vector2d(cache.dff[k][i].data());
Eigen::MatrixXd fk;
fct.eval(*this, u, k, fk, fake);
- double detJ = mem.getDetJ(k);
- Eigen::Matrix2d const &J = mem.getJinv(k);
// Elementary stiffness matrix
- for (size_t i = 0; i < nodes.size(); ++i)
- {
- std::vector<double> &dffi = cache.dff[k][i];
- factor0 = J * Eigen::Vector2d(dffi.data());
- factor1 = fk * factor0;
-
- for (size_t j = 0; j < nodes.size(); ++j)
- {
- std::vector<double> &dffj = cache.dff[k][j];
- factor2 = J * Eigen::Vector2d(dffj.data());
- K(i, j) += factor1.dot(factor2) * gauss.w[k] * detJ;
- }
- }
+ K += (fk * J * dN).transpose() * J * dN * gauss.w[k] * mem.getDetJ(k);
}
}
}
@@ -208,35 +193,33 @@ void Tri3::buildK(Eigen::MatrixXd &K, std::vector<double> const &u, Fct2 const &
* K_ij(6,6) = int eps_k H_ijkl eps_l dV = B(6,3)H(3,3)B(3,6)
* @authors Adrien Crovato
*/
-void Tri3::buildK_linElastic(Eigen::MatrixXd &K, Eigen::MatrixXd const &H)
+void Tri3::buildK_mechanical(Eigen::MatrixXd &K, Eigen::MatrixXd const &H)
{
CacheTri3 &cache = getCache();
GaussTri3 &gauss = cache.gauss;
MemTri3 &mem = getMem();
+ // Gauss integration
K = Eigen::MatrixXd::Zero(6, 6);
- Eigen::Vector2d factori;
for (size_t k = 0; k < gauss.p.size(); ++k)
{
// Jacobian inverse
- double detJ = mem.getDetJ(k);
Eigen::Matrix2d const &invJ = mem.getJinv(k);
-
- // Fill B matrix
+ // Fill B matrix (shape functions)
Eigen::MatrixXd B = Eigen::MatrixXd::Zero(3, 6);
for (size_t i = 0; i < nodes.size(); ++i)
{
// Compute [Ni,x Ni,y]
std::vector<double> &dffi = cache.dff[k][i];
- factori = invJ * Eigen::Vector2d(dffi.data());
- B(0, 2 * i) = factori(0);
- B(1, 2 * i + 1) = factori(1);
- B(2, 2 * i + 1) = factori(0);
- B(2, 2 * i) = factori(1);
+ Eigen::Vector2d dN = invJ * Eigen::Vector2d(dffi.data());
+ B(0, 2 * i) = dN(0);
+ B(1, 2 * i + 1) = dN(1);
+ B(2, 2 * i + 1) = dN(0);
+ B(2, 2 * i) = dN(1);
}
// Compute stiffness matrix
- K += B.transpose() * H * B * detJ * gauss.w[k];
+ K += B.transpose() * H * B * gauss.w[k] * mem.getDetJ(k);
}
}
@@ -252,30 +235,20 @@ void Tri3::buildDK(Eigen::MatrixXd &K, std::vector<double> const &u, Fct0 const
GaussTri3 &gauss = cache.gauss;
MemTri3 &mem = getMem();
+ // Gauss integration
K = Eigen::Matrix3d::Zero();
- Eigen::Vector2d factor1;
- Eigen::Vector2d factor2;
for (size_t k = 0; k < gauss.p.size(); ++k)
{
- // "derivative" of f, gradient and inverse Jacobian
- double dfk = fct.evalD(*this, u, k);
+ // Jacobian inverse, shape functions and gradient
+ Eigen::Matrix2d const &J = mem.getJinv(k);
+ Eigen::Matrix<double, 2, 3> dN;
+ for (size_t i = 0; i < nodes.size(); ++i)
+ dN.col(i) = Eigen::Vector2d(cache.dff[k][i].data());
Eigen::VectorXd gradu;
this->computeGrad(u, gradu, k);
- double detJ = mem.getDetJ(k);
- Eigen::Matrix2d const &J = mem.getJinv(k);
- // computation of K
- for (size_t i = 0; i < nodes.size(); ++i)
- {
- std::vector<double> &dffi = cache.dff[k][i];
- factor1 = J * Eigen::Vector2d(dffi.data());
- for (size_t j = 0; j < nodes.size(); ++j)
- {
- std::vector<double> &dffj = cache.dff[k][j];
- factor2 = J * Eigen::Vector2d(dffj.data());
- K(i, j) += gradu.dot(factor1) * gradu.dot(factor2) * gauss.w[k] * detJ * dfk;
- }
- }
+ // Elementary stiffness matrix
+ K += (gradu.transpose() * J * dN).transpose() * gradu.transpose() * J * dN * fct.evalD(*this, u, k) * gauss.w[k] * mem.getDetJ(k);
}
}
@@ -291,13 +264,10 @@ void Tri3::builds(Eigen::VectorXd &s, std::vector<double> const &u, Fct0 const &
GaussTri3 &gauss = cache.gauss;
MemTri3 &mem = getMem();
+ // Gauss integration
s = Eigen::Vector3d::Zero();
for (size_t k = 0; k < gauss.p.size(); ++k)
- {
- double sk = f.eval(*this, u, k);
- double detJ = mem.getDetJ(k);
- s += Eigen::Vector3d(cache.ff[k].data()) * sk * gauss.w[k] * detJ; // s = s + sk.wG.detJ
- }
+ s += Eigen::Vector3d(cache.ff[k].data()) * f.eval(*this, u, k) * gauss.w[k] * mem.getDetJ(k);
}
/**
@@ -312,13 +282,11 @@ void Tri3::builds(Eigen::VectorXd &s, std::vector<double> const &u, Fct1 const &
GaussTri3 &gauss = cache.gauss;
MemTri3 &mem = getMem();
+ // Gauss integration
s = Eigen::Vector3d::Zero();
for (size_t k = 0; k < gauss.p.size(); ++k)
- {
- double sk = f.eval(*this, u, k) * normal();
- double detJ = mem.getDetJ(k);
- s += Eigen::Vector3d(cache.ff[k].data()) * sk * gauss.w[k] * detJ; // s = s + sk.wG.detJ
- }
+ s += Eigen::Vector3d(cache.ff[k].data()) * (f.eval(*this, u, k) * normal()) * gauss.w[k] * mem.getDetJ(k);
+
}
/**
@@ -333,30 +301,26 @@ void Tri3::buildDs(Eigen::VectorXd &s, std::vector<double> const &u, Fct0 const
GaussTri3 &gauss = cache.gauss;
MemTri3 &mem = getMem();
+ // Gauss integration
s = Eigen::Vector3d::Zero();
- Eigen::Vector2d factor1;
for (size_t k = 0; k < gauss.p.size(); ++k)
{
+ // Shape functions and gradient
+ Eigen::Matrix<double, 2, 3> dN;
+ for (size_t i = 0; i < nodes.size(); ++i)
+ dN.col(i) = Eigen::Vector2d(cache.dff[k][i].data());
Eigen::VectorXd gradu;
this->computeGrad(u, gradu, k);
- double fk = f.eval(*this, u, k);
- double detJ = mem.getDetJ(k);
- Eigen::Matrix2d const &J = mem.getJinv(k);
- // s(i) = f * gradu * gradN(i)
- for (size_t i = 0; i < nodes.size(); ++i)
- {
- std::vector<double> &dffi = cache.dff[k][i];
- factor1 = J * Eigen::Vector2d(dffi.data());
- s(i) += gradu.dot(factor1) * gauss.w[k] * detJ * fk;
- }
+ // Elementary flux vector
+ s += gradu.transpose() * mem.getJinv(k) * dN * f.eval(*this, u, k) * gauss.w[k] * mem.getDetJ(k);
}
}
/**
- * @brief Compute something
+ * @brief Compute mass matrix
*
- * S_ij(3,3) = int N_i*N_j dV (from Quad4)
+ * S_ij(3,3) = int N_i*N_j dV
* @authors Romain Boman
*/
void Tri3::buildS(Eigen::MatrixXd &S)
@@ -365,6 +329,7 @@ void Tri3::buildS(Eigen::MatrixXd &S)
GaussTri3 &gauss = cache.gauss;
MemTri3 &mem = getMem();
+ // Gauss integration
S = Eigen::Matrix3d::Zero();
for (size_t k = 0; k < gauss.p.size(); ++k)
{
@@ -372,10 +337,14 @@ void Tri3::buildS(Eigen::MatrixXd &S)
double detJ = mem.getDetJ(k);
std::vector<double> &ff = cache.ff[k];
- S += Eigen::Vector3d(ff.data()) * Eigen::RowVector3d(ff.data()) * gauss.w[k] * detJ; // S = S + {ff}^T [ff] * wG * detJ
+ S += Eigen::Vector3d(ff.data()) * Eigen::RowVector3d(ff.data()) * gauss.w[k] * detJ;
}
}
+/**
+ * @brief API for build methods
+ * @authors Romain Boman
+ */
void Tri3::build(MATTYPE type, Eigen::MatrixXd &A, std::vector<double> const &u, Fct0 const &fct)
{
switch (type)
@@ -402,7 +371,7 @@ void Tri3::computeGrad(std::vector<double> const &u, Eigen::VectorXd &grad, size
CacheTri3 &cache = getCache();
MemTri3 &mem = getMem();
- // local matrix of sf derivatives, solution vector and Jacobian
+ // local matrix of sf derivatives and solution vector
Eigen::MatrixXd gradN(2, nodes.size());
Eigen::Vector3d ue;
for (size_t i = 0; i < nodes.size(); ++i)
@@ -410,9 +379,8 @@ void Tri3::computeGrad(std::vector<double> const &u, Eigen::VectorXd &grad, size
gradN.col(i) = Eigen::Vector2d(cache.dff[k][i].data());
ue(i) = u[nodes[i]->row];
}
- Eigen::Matrix2d const &Jinv = mem.getJinv(k);
// Gradient in global axis: inv(J)_ij * d_jN_i * ue_i
- grad = Jinv * gradN * ue;
+ grad = mem.getJinv(k) * gradN * ue;
}
/**
@@ -427,7 +395,10 @@ double Tri3::computeGrad2(std::vector<double> const &u, size_t k)
}
/**
- * calculation of int_V -gradu.fct dV
+ * @brief Compute flux
+ *
+ * qV_i(2) = int -gradu * f dV
+ * @authors Romain Boman
*/
void Tri3::computeqV(std::vector<double> const &u, Eigen::VectorXd &qV, Fct0 const &fct)
{
@@ -435,21 +406,24 @@ void Tri3::computeqV(std::vector<double> const &u, Eigen::VectorXd &qV, Fct0 con
GaussTri3 &gauss = cache.gauss;
MemTri3 &mem = getMem();
+ // Gauss integration
qV = Eigen::Vector2d::Zero();
for (size_t k = 0; k < gauss.p.size(); ++k)
{
- // Gradient function and Jacobian evaluation
+ // Gradient
Eigen::VectorXd gradk;
computeGrad(u, gradk, k);
- double fk = fct.eval(*this, u, k);
- double detJ = mem.getDetJ(k);
- qV -= gradk * fk * gauss.w[k] * detJ; // qV = qV + fk.wG.detJ
+ // Elementary flux
+ qV -= gradk * fct.eval(*this, u, k) * gauss.w[k] * mem.getDetJ(k);
}
}
/**
- * calculation of int_V flux grad N_j dV
+ * @brief Compute flux
+ *
+ * qV_i(3) = int [f] gradu N_i dV
+ * @authors Romain Boman
*/
void Tri3::computeqint(std::vector<double> const &u, Eigen::VectorXd &qV, Fct2 const &fct)
{
@@ -457,32 +431,29 @@ void Tri3::computeqint(std::vector<double> const &u, Eigen::VectorXd &qV, Fct2 c
GaussTri3 &gauss = cache.gauss;
MemTri3 &mem = getMem();
+ // Gauss integration
qV = Eigen::Vector3d::Zero();
for (size_t k = 0; k < gauss.p.size(); ++k)
{
- // Gradient and flux evaluation
+ // Shape functions, gradient and flux evaluation
+ Eigen::Matrix<double, 2, 3> dN;
+ for (size_t i = 0; i < nodes.size(); ++i)
+ dN.col(i) = Eigen::Vector2d(cache.dff[k][i].data());
Eigen::VectorXd gradk;
computeGrad(u, gradk, k);
Eigen::MatrixXd fk(2, 2);
fct.eval(*this, u, k, fk, false);
- double detJ = mem.getDetJ(k);
- Eigen::Vector2d flux;
- flux = fk * gradk;
- // calcul du gradient des fcts de forme et produit scalaire
- Eigen::Vector2d factor1;
- Eigen::Matrix2d const &J = mem.getJinv(k);
- for (size_t i = 0; i < nodes.size(); ++i)
- {
- std::vector<double> &dffi = cache.dff[k][i];
- factor1 = J * Eigen::Vector2d(dffi.data());
- qV(i) += factor1.dot(flux) * gauss.w[k] * detJ;
- }
+ // Elementary flux vector
+ qV += (fk *gradk).transpose() * mem.getJinv(k) * dN * gauss.w[k] * mem.getDetJ(k);
}
}
/**
- * calculation of int_V -[fct2]gradu dV
+ * @brief Compute flux
+ *
+ * qV_i(3) = int -[f] gradu dV
+ * @authors Romain Boman
*/
void Tri3::computeqV(std::vector<double> const &u, Eigen::VectorXd &qV, Fct2 const &fct)
{
@@ -490,6 +461,7 @@ void Tri3::computeqV(std::vector<double> const &u, Eigen::VectorXd &qV, Fct2 con
GaussTri3 &gauss = cache.gauss;
MemTri3 &mem = getMem();
+ // Gauss integration
qV = Eigen::Vector2d::Zero();
for (size_t k = 0; k < gauss.p.size(); ++k)
{
@@ -498,16 +470,15 @@ void Tri3::computeqV(std::vector<double> const &u, Eigen::VectorXd &qV, Fct2 con
computeGrad(u, gradk, k);
Eigen::MatrixXd fk(2, 2);
fct.eval(*this, u, k, fk, false);
- double detJ = mem.getDetJ(k);
- qV -= fk * gradk * gauss.w[k] * detJ; // qV = qV + fk.wG.detJ
+ qV -= fk * gradk * gauss.w[k] * mem.getDetJ(k);
}
}
/**
- * @brief Compute intergal of function
+ * @brief Compute intergal of scalar function
*
- * I = int fct dV
+ * I = int f dV
* @authors Romain Boman
*/
double Tri3::computeV(std::vector<double> const &u, Fct0 const &fct)
@@ -516,34 +487,42 @@ double Tri3::computeV(std::vector<double> const &u, Fct0 const &fct)
GaussTri3 &gauss = cache.gauss;
MemTri3 &mem = getMem();
+ // Gauss integration
double V = 0.0;
for (size_t k = 0; k < gauss.p.size(); ++k)
- {
- // Function evaluation
- double fk = fct.eval(*this, u, k);
- double detJ = mem.getDetJ(k);
- V += fk * gauss.w[k] * detJ;
- }
+ V += fct.eval(*this, u, k) * gauss.w[k] * mem.getDetJ(k);
return V;
}
+
+/**
+ * @brief Compute intergal of matrix function
+ *
+ * I = int [f] dV
+ * @authors Romain Boman
+ */
void Tri3::computeV(std::vector<double> const &u, Eigen::MatrixXd &out, Fct2 const &fct)
{
CacheTri3 &cache = getCache();
GaussTri3 &gauss = cache.gauss;
MemTri3 &mem = getMem();
- out = Eigen::MatrixXd::Zero(2, 2); // todo out should be resized to match fk size
+ // Gauss integration
+ out = Eigen::MatrixXd::Zero(2, 2); /* @todo out should be resized to match fk size */
for (size_t k = 0; k < gauss.p.size(); ++k)
{
// Function evaluation
Eigen::MatrixXd fk;
fct.eval(*this, u, k, fk, false);
- double detJ = mem.getDetJ(k);
- out += fk * gauss.w[k] * detJ; // [out] = [out] + [fk].wG.detJ
+ out += fk * gauss.w[k] * mem.getDetJ(k);
}
}
+/**
+ * @brief ???
+ * @todo never used, should be removed
+ * @authors Romain Boman
+ */
void Tri3::iso_coordinates(Pt x, double eta, double dzeta)
{
Eigen::Vector2d iso, b;
@@ -602,11 +581,10 @@ void Tri3::iso_coordinates(Pt x, double eta, double dzeta)
}
/**
- * @brief Return the element volume
+ * @brief Return the element surface
* @authors Romain Boman
*/
double Tri3::computeV()
{
- MemTri3 &mem = getMem();
- return mem.getVol();
+ return getMem().getVol();
}
diff --git a/tbox/src/wTri3.h b/tbox/src/wTri3.h
index d6b843c1..2102c993 100644
--- a/tbox/src/wTri3.h
+++ b/tbox/src/wTri3.h
@@ -40,7 +40,7 @@ class TBOX_API Tri3 : public Element
virtual void build(MATTYPE type, Eigen::MatrixXd &A, std::vector<double> const &u, Fct0 const &fct);
virtual void buildK(Eigen::MatrixXd &K, std::vector<double> const &u, Fct0 const &fct) override;
virtual void buildK(Eigen::MatrixXd &K, std::vector<double> const &u, Fct2 const &fct, bool fake) override;
- virtual void buildK_linElastic(Eigen::MatrixXd &K, Eigen::MatrixXd const &H) override;
+ virtual void buildK_mechanical(Eigen::MatrixXd &K, Eigen::MatrixXd const &H) override;
virtual void builds(Eigen::VectorXd &s, std::vector<double> const &u, Fct0 const &f) override;
virtual void builds(Eigen::VectorXd &s, std::vector<double> const &u, Fct1 const &f) override;
virtual void buildDK(Eigen::MatrixXd &K, std::vector<double> const &u, Fct0 const &fct) override;
--
GitLab