diff --git a/flow/src/wF0El.cpp b/flow/src/wF0El.cpp
index 45a27a4c9a9cf0bf07eb4a03522be0241d7896ec..5643219538e53a449da43f1c0d6ca8f020c465dc 100644
--- a/flow/src/wF0El.cpp
+++ b/flow/src/wF0El.cpp
@@ -24,7 +24,7 @@ using namespace flow;
*/
double F0ElRho::eval(Element const &e, std::vector<double> const &u, size_t k) const
{
- double gradU = e.computeGrad(u, k).norm(); // norm of velocity
+ double gradU = e.computeGradient(u, k).norm(); // norm of velocity
if (gradU < gradUC) // true density
{
@@ -45,7 +45,7 @@ double F0ElRho::eval(Element const &e, std::vector<double> const &u, size_t k) c
*/
double F0ElRho::evalD(Element const &e, std::vector<double> const &u, size_t k) const
{
- double gradU = e.computeGrad(u, k).norm(); // norm of velocity
+ double gradU = e.computeGradient(u, k).norm(); // norm of velocity
if (gradU < gradUC) // true density
{
@@ -92,7 +92,8 @@ void F0ElRhoL::write(std::ostream &out) const
*/
double F0ElMach::eval(Element const &e, std::vector<double> const &u, size_t k) const
{
- double gradU2 = e.computeGrad2(u, k); // velocity squared
+ Eigen::VectorXd gradU = e.computeGradient(u, k); // velocity
+ double gradU2 = gradU.squaredNorm(); // velocity squared
double a2 = 1 / (mInf * mInf) + 0.5 * (gamma - 1) * (1 - gradU2); // speed of sound squared
return sqrt(gradU2) / sqrt(a2);
@@ -103,7 +104,8 @@ double F0ElMach::eval(Element const &e, std::vector<double> const &u, size_t k)
*/
double F0ElMach::evalD(Element const &e, std::vector<double> const &u, size_t k) const
{
- double gradU2 = e.computeGrad2(u, k);
+ Eigen::VectorXd gradU = e.computeGradient(u, k); // velocity
+ double gradU2 = gradU.squaredNorm(); // velocity squared
double a2 = 1 / (mInf * mInf) + 0.5 * (gamma - 1) * (1 - gradU2);
return 1 / sqrt(gradU2 * a2) + 0.5 * (gamma - 1) * sqrt(gradU2) / sqrt(a2 * a2 * a2); // has to be multiplied by grad u to recover true derivative
@@ -140,7 +142,8 @@ void F0ElMachL::write(std::ostream &out) const
*/
double F0ElCp::eval(Element const &e, std::vector<double> const &u, size_t k) const
{
- double gradU2 = e.computeGrad2(u, k); // velocity squared
+ Eigen::VectorXd gradU = e.computeGradient(u, k); // velocity
+ double gradU2 = gradU.squaredNorm(); // velocity squared
//particularized: 2 / (gamma * mInf * mInf) * (pow(1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradU2), gamma / (gamma - 1)) - 1);
double a = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradU2);
@@ -155,7 +158,8 @@ double F0ElCp::eval(Element const &e, std::vector<double> const &u, size_t k) co
*/
double F0ElCp::evalD(Element const &e, std::vector<double> const &u, size_t k) const
{
- double gradU2 = e.computeGrad2(u, k); // velocity squared
+ Eigen::VectorXd gradU = e.computeGradient(u, k); // velocity
+ double gradU2 = gradU.squaredNorm(); // velocity squared
//particularized: -2 * pow(1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradU2), 1 / (gamma - 1));
double a = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradU2);
@@ -175,7 +179,9 @@ void F0ElCp::write(std::ostream &out) const
*/
double F0ElCpL::eval(Element const &e, std::vector<double> const &u, size_t k) const
{
- return 1 - e.computeGrad2(u, k);
+ Eigen::VectorXd gradU = e.computeGradient(u, k); // velocity
+ double gradU2 = gradU.squaredNorm(); // velocity squared
+ return 1 - gradU2;
}
double F0ElCpL::evalD(Element const &e, std::vector<double> const &u, size_t k) const
{
diff --git a/flow/src/wKuttaResidual.cpp b/flow/src/wKuttaResidual.cpp
index 59c2010a5453f0d76a357af6d4f6b32e23aa9b57..d7adcac995d1b079033d0e95dcab1272eaae953d 100644
--- a/flow/src/wKuttaResidual.cpp
+++ b/flow/src/wKuttaResidual.cpp
@@ -40,7 +40,7 @@ Eigen::MatrixXd KuttaResidual::buildFixed(KuttaElement const &ke, std::vector<do
size_t nRow = ke.nRow;
// Velocity
- Eigen::RowVectorXd Aup = ke.volE->computeGrad(phi, 0).transpose() * ke.volE->getVMem().getJinv(0) * ke.volE->getVCache().getDsf(0);
+ Eigen::RowVectorXd Aup = ke.volE->computeGradient(phi, 0).transpose() * ke.volE->getVMem().getJinv(0) * ke.volE->getVCache().getDsf(0);
// Build
Eigen::MatrixXd A = Eigen::MatrixXd::Zero(nRow, ke.nCol);
@@ -72,7 +72,8 @@ Eigen::VectorXd KuttaResidual::build(KuttaElement const &ke, std::vector<double>
size_t nRow = ke.nRow;
// Velocity
- double dPhi2 = ke.volE->computeGrad2(phi, 0);
+ Eigen::VectorXd dPhi = ke.volE->computeGradient(phi, 0);
+ double dPhi2 = dPhi.squaredNorm();
// Build
Eigen::VectorXd b = Eigen::VectorXd::Zero(nRow);
@@ -125,7 +126,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> KuttaResidual::buildGradientMesh(Ku
size_t nSur = ke.surE->nodes.size();
// velocity
- Eigen::VectorXd dPhi = ke.volE->computeGrad(phi, 0);
+ Eigen::VectorXd dPhi = ke.volE->computeGradient(phi, 0);
// Build
Eigen::MatrixXd Av = Eigen::MatrixXd::Zero(nRow, nDim * nCol);
diff --git a/flow/src/wLoadFunctional.cpp b/flow/src/wLoadFunctional.cpp
index a251f54f689d3fe9ec25a4f6112e467844faf324..feefdaf9b2fdfc89b04b9fcb2dec2dd8a58acd17 100644
--- a/flow/src/wLoadFunctional.cpp
+++ b/flow/src/wLoadFunctional.cpp
@@ -65,7 +65,7 @@ Eigen::MatrixXd LoadFunctional::buildGradientFlow(Element const &e, Element cons
Mem &mem = e.getVMem();
// Gradient of pressure force coefficient in associated volume
- Eigen::MatrixXd dfcp = e.normal() * cp.evalD(eV, phi, 0) * eV.computeGrad(phi, 0).transpose() * eV.getVMem().getJinv(0) * eV.getVCache().getDsf(0);
+ Eigen::MatrixXd dfcp = e.normal() * cp.evalD(eV, phi, 0) * eV.computeGradient(phi, 0).transpose() * eV.getVMem().getJinv(0) * eV.getVCache().getDsf(0);
// Integrate cp on the surface
Eigen::MatrixXd A = Eigen::MatrixXd::Zero(3 * e.nodes.size(), eV.nodes.size());
@@ -99,7 +99,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> LoadFunctional::buildGradientMesh(E
Eigen::Vector3d nrm = e.normal(); // unit normal
// Gradients
double dCp = cp.evalD(eV, phi, 0); // pressure coefficient
- Eigen::VectorXd dPhi = eV.computeGrad(phi, 0); // potential
+ 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
diff --git a/flow/src/wNewton.cpp b/flow/src/wNewton.cpp
index 4ac08c51708639d0c16c8e06b8a2aabda7e75665..71a50f75ca6e842d13d9e4d35807aa3d217c4520 100644
--- a/flow/src/wNewton.cpp
+++ b/flow/src/wNewton.cpp
@@ -502,7 +502,7 @@ void Newton::findUpwd()
{
tbb::parallel_do(pbl->medium->adjMap.begin(), pbl->medium->adjMap.end(), [&](std::pair<Element *, std::vector<Element *>> &p) {
// Opposite velocity vector
- Eigen::VectorXd vel_ = -p.first->computeGrad(phi, 0);
+ Eigen::VectorXd vel_ = -p.first->computeGradient(phi, 0);
vel_.normalize();
// Best alignement strategy
// -> find (cgU-cgE) \dot vel_ closest to one
diff --git a/flow/src/wPotentialResidual.cpp b/flow/src/wPotentialResidual.cpp
index ff2054788bdefe5fd12d4ebcfcbdaffb128f1552..f98adbc55f97dd7e3126fb12c4e7c300239e2570 100644
--- a/flow/src/wPotentialResidual.cpp
+++ b/flow/src/wPotentialResidual.cpp
@@ -65,7 +65,7 @@ Eigen::VectorXd PotentialResidual::build(Element const &e, Element const &eU, st
// 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.computeGrad(phi, k) * gauss.getW(k) * mem.getDetJ(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);
// Supersonic contribution
double mach = fluid.mach.eval(e, phi, 0);
@@ -77,7 +77,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.computeGrad(phi, k) * gauss.getW(k) * mem.getDetJ(k);
+ b += mu * rhoU * (mem.getJinv(k) * cache.getDsf(k)).transpose() * e.computeGradient(phi, k) * gauss.getW(k) * mem.getDetJ(k);
}
return b;
}
@@ -107,7 +107,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientFlo
double wdj = gauss.getW(k) * mem.getDetJ(k);
// Shape functions and solution gradient
Eigen::MatrixXd const &dSf = mem.getJinv(k) * cache.getDsf(k);
- Eigen::VectorXd dPhi = e.computeGrad(phi, k);
+ Eigen::VectorXd dPhi = e.computeGradient(phi, k);
// rho * grad_phi*grad_psi + drho * grad_phi*grad_psi
A1 += fluid.rho.eval(e, phi, k) * dSf.transpose() * dSf * wdj;
@@ -123,7 +123,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientFlo
double mu = muC * (1 - (mCO * mCO) / (mach * mach));
double dmu = 2 * muC * mCO * mCO / (mach * mach * mach);
// density and gradients on upwind elements
- Eigen::VectorXd dPhiU = eU.computeGrad(phi, 0);
+ 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);
@@ -137,7 +137,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientFlo
double wdj = gauss.getW(k) * mem.getDetJ(k);
// Shape functions and solution gradient
Eigen::MatrixXd const &dSf = mem.getJinv(k) * cache.getDsf(k);
- Eigen::VectorXd dPhi = e.computeGrad(phi, k);
+ Eigen::VectorXd dPhi = e.computeGradient(phi, k);
// mu*drhoU*grad_phi*grad_psi
A2 += mu * dRhoU * dSf.transpose() * dPhi * dPhiU.transpose() * dSfU * wdj;
@@ -173,7 +173,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientMes
for (size_t k = 0; k < gauss.getN(); ++k)
{
// Density and gradients
- Eigen::VectorXd dPhi = e.computeGrad(phi, k);
+ Eigen::VectorXd dPhi = e.computeGradient(phi, k);
double rho = fluid.rho.eval(e, phi, k);
double dRho = fluid.rho.evalD(e, phi, k);
Eigen::MatrixXd const &dSf = cache.getDsf(k);
@@ -205,7 +205,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientMes
double mu = muC * (1 - (mCO * mCO) / (mach * mach));
double dmu = 2 * muC * mCO * mCO / (mach * mach * mach);
// upwind density and gradient
- Eigen::VectorXd dPhiU = eU.computeGrad(phi, 0);
+ 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);
@@ -217,7 +217,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientMes
for (size_t k = 0; k < gauss.getN(); ++k)
{
// Density, Mach and Gradients
- Eigen::VectorXd dPhi = e.computeGrad(phi, k);
+ Eigen::VectorXd dPhi = e.computeGradient(phi, k);
double rho = fluid.rho.eval(e, phi, k);
double dM = fluid.mach.evalD(e, phi, k);
Eigen::MatrixXd const &dSf = cache.getDsf(k);
diff --git a/flow/src/wSolver.cpp b/flow/src/wSolver.cpp
index c90dda9f650ce5df742dbd6e46e7b580b62e7aa1..d44e29dcfdaa7af01618a2e658dfe728fec40961 100644
--- a/flow/src/wSolver.cpp
+++ b/flow/src/wSolver.cpp
@@ -228,7 +228,7 @@ void Solver::computeFlow()
U[n->row] = Eigen::Vector3d::Zero();
for (auto e : p.second)
{
- Eigen::VectorXd grad = e->computeGrad(phi, 0);
+ Eigen::VectorXd grad = e->computeGradient(phi, 0);
for (size_t i = 0; i < grad.size(); ++i)
U[n->row](i) += grad(i);
}
diff --git a/flow/src/wWakeResidual.cpp b/flow/src/wWakeResidual.cpp
index 934b69a295aa60cfb05d0e9223c6c1541db0bdb6..0b111e233d2e82a703972e027aa134bbe739a014 100644
--- a/flow/src/wWakeResidual.cpp
+++ b/flow/src/wWakeResidual.cpp
@@ -55,8 +55,8 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildFixed(WakeElemen
Eigen::VectorXd dSfp(nRow);
dSfp = 0.5 * sqrt(sMem.getVol()) * dSf * vUpJ.transpose() * vi;
// Velocity
- Eigen::RowVectorXd Aup = we.volUpE->computeGrad(phi, 0).transpose() * vUpJ * vUpDsf;
- Eigen::RowVectorXd Alw = we.volLwE->computeGrad(phi, 0).transpose() * vLwJ * vLwDsf;
+ Eigen::RowVectorXd Aup = we.volUpE->computeGradient(phi, 0).transpose() * vUpJ * vUpDsf;
+ Eigen::RowVectorXd Alw = we.volLwE->computeGradient(phi, 0).transpose() * vLwJ * vLwDsf;
// Build
Eigen::MatrixXd Au = Eigen::MatrixXd::Zero(nRow, we.nColUp);
@@ -101,8 +101,10 @@ std::tuple<Eigen::VectorXd, Eigen::VectorXd> WakeResidual::build(WakeElement con
Eigen::VectorXd dSfp(nRow);
dSfp = 0.5 * sqrt(sMem.getVol()) * dSf * we.volUpE->getVMem().getJinv(0).transpose() * vi;
// Velocity squared
- double dPhi2u = we.volUpE->computeGrad2(phi, 0);
- double dPhi2l = we.volLwE->computeGrad2(phi, 0);
+ Eigen::VectorXd dPhiU = we.volUpE->computeGradient(phi, 0);
+ Eigen::VectorXd dPhiL = we.volLwE->computeGradient(phi, 0);
+ double dPhi2u = dPhiU.squaredNorm();
+ double dPhi2l = dPhiL.squaredNorm();
// Build
Eigen::VectorXd bu = Eigen::VectorXd::Zero(nRow);
@@ -193,8 +195,8 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buil
}
}
// Gradients of potential
- Eigen::VectorXd dPhiUp = we.volUpE->computeGrad(phi, 0);
- Eigen::VectorXd dPhiLw = we.volLwE->computeGrad(phi, 0);
+ Eigen::VectorXd dPhiUp = we.volUpE->computeGradient(phi, 0);
+ Eigen::VectorXd dPhiLw = we.volLwE->computeGradient(phi, 0);
// Build
Eigen::MatrixXd Au = Eigen::MatrixXd::Zero(nRow, nDim * nColUp);
diff --git a/heat/src/heat.h b/heat/src/heat.h
index 479b9945ae3087528f31e468b78f6a9390cb5296..5084a8f82a0c8da88c216c62e8b4ff5a3ecfc785 100644
--- a/heat/src/heat.h
+++ b/heat/src/heat.h
@@ -48,6 +48,8 @@ class Extractor;
class NodePair;
+class HeatTerm;
+
// not used yet
class DisplayHook;
}; // namespace heat
diff --git a/heat/src/wHeatTerm.cpp b/heat/src/wHeatTerm.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..057e139c461a5f2b5d2b2a240d158a8ef448db9b
--- /dev/null
+++ b/heat/src/wHeatTerm.cpp
@@ -0,0 +1,154 @@
+/*
+ * 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 "wHeatTerm.h"
+
+#include "wElement.h"
+#include "wCache.h"
+#include "wGauss.h"
+#include "wMem.h"
+#include "wFct0.h"
+#include "wFct2.h"
+
+using namespace heat;
+using namespace tbox;
+
+/**
+ * @brief Build volume terms for the heat equation on one element
+ */
+Eigen::MatrixXd HeatTerm::build(Element const &e, std::vector<double> const &u, Fct2 const &f, bool fake)
+{
+ // Get precomputed values
+ Cache &cache = e.getVCache();
+ Gauss &gauss = cache.getVGauss();
+
+ Eigen::MatrixXd K = Eigen::MatrixXd::Zero(e.nodes.size(), e.nodes.size());
+ if (fake) // fake run - fill MPI job list
+ {
+ for (size_t k = 0; k < gauss.getN(); ++k)
+ {
+ Eigen::MatrixXd fk;
+ f.eval(e, u, k, fk, fake);
+ }
+ }
+ 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 &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);
+ }
+ }
+ return K;
+}
+
+/**
+ * @brief Build surface/volume flux on one element
+ */
+Eigen::VectorXd HeatTerm::build(Element const &e, std::vector<double> const &u, Fct0 const &f)
+{
+ // 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);
+ return s;
+}
+
+/**
+ * @brief Build volume flux of the heat equation on one element
+ * @todo not used
+ */
+Eigen::VectorXd HeatTerm::build2(Element const &e, std::vector<double> const &u, Fct2 const &f)
+{
+ // 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());
+ for (size_t k = 0; k < gauss.getN(); ++k)
+ {
+ // Shape functions, gradient and flux evaluation
+ Eigen::MatrixXd const &dff = cache.getDsf(k);
+ Eigen::MatrixXd fk(2, 2);
+ 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);
+ }
+ return qV;
+}
+
+/**
+ * @brief Compute (integrate) volume flux
+ */
+Eigen::VectorXd HeatTerm::computeFlux(Element const &e, std::vector<double> const &u, Fct2 const &f)
+{
+ // 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
+ for (size_t k = 0; k < gauss.getN(); ++k)
+ {
+ // Gradient and matrix function evaluation
+ Eigen::VectorXd gradk = e.computeGradient(u, k);
+ Eigen::MatrixXd fk(2, 2);
+ f.eval(e, u, k, fk, false);
+
+ qV -= fk * gradk * gauss.getW(k) * mem.getDetJ(k);
+ }
+ return qV;
+}
+
+/**
+ * @brief Compute (integrate) matrix
+ */
+Eigen::MatrixXd HeatTerm::computeMatrix(Element const &e, std::vector<double> const &u, Fct2 const &f)
+{
+ // 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
+ for (size_t k = 0; k < gauss.getN(); ++k)
+ {
+ // Function evaluation
+ Eigen::MatrixXd fk;
+ f.eval(e, u, k, fk, false);
+
+ out += fk * gauss.getW(k) * mem.getDetJ(k);
+ }
+ return out;
+}
diff --git a/heat/src/wHeatTerm.h b/heat/src/wHeatTerm.h
new file mode 100644
index 0000000000000000000000000000000000000000..3bbfada2079b4f9522c599c208f0baca3b57a5fc
--- /dev/null
+++ b/heat/src/wHeatTerm.h
@@ -0,0 +1,45 @@
+/*
+ * 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 WHEATTERM_H
+#define WHEATTERM_H
+
+#include "heat.h"
+
+#include <vector>
+#include <Eigen/Dense>
+
+namespace heat
+{
+
+/**
+ * @brief Formulation of heat terms
+ */
+class HEAT_API HeatTerm
+{
+public:
+ // Volume
+ static Eigen::MatrixXd build(tbox::Element const &e, std::vector<double> const &u, tbox::Fct2 const &f, bool fake);
+ static Eigen::VectorXd build2(tbox::Element const &e, std::vector<double> const &u, tbox::Fct2 const &f);
+ // Boundary and Source
+ static Eigen::VectorXd build(tbox::Element const &e, std::vector<double> const &u, tbox::Fct0 const &f);
+ // Functionals
+ static Eigen::VectorXd computeFlux(tbox::Element const &e, std::vector<double> const &u, tbox::Fct2 const &f);
+ static Eigen::MatrixXd computeMatrix(tbox::Element const &e, std::vector<double> const &u, tbox::Fct2 const &f);
+};
+
+} // namespace heat
+#endif //WHEATTERM_H
diff --git a/heat/src/wSolver.cpp b/heat/src/wSolver.cpp
index 8600d1ea6c25e5b0c224136cfa68501c53228027..685cabd12719e02da28920275a06d8cc604eee26 100644
--- a/heat/src/wSolver.cpp
+++ b/heat/src/wSolver.cpp
@@ -29,6 +29,7 @@
#include "wDirichlet.h"
#include "wDisplayHook.h"
#include "wPeriodic.h"
+#include "wHeatTerm.h"
#include "wResults.h"
#include "wMshExport.h"
#include "wLinearSolver.h"
@@ -263,14 +264,14 @@ void Solver::start(MshExport *mshWriter)
//std::cout << "processing element #" << e->no << "\n";
// (flux moyen . volume) sur l'element
- Eigen::VectorXd qV = e->computeqV(T1, mat->k);
+ Eigen::VectorXd qV = HeatTerm::computeFlux(*e, T1, mat->k);
double V = e->computeV();
if (save)
{
// gradient of T at GP 0
- Eigen::VectorXd grad = e->computeGrad(T1, 0);
+ Eigen::VectorXd grad = e->computeGradient(T1, 0);
int i = e->no - 1;
gradT_x[i] = grad(0);
gradT_y[i] = grad(1);
@@ -280,7 +281,7 @@ void Solver::start(MshExport *mshWriter)
kgradT[i] = Eigen::Vector3d(qV(0) / V, qV(1) / V, 0);
// mean k over the element
- Eigen::MatrixXd kmean = e->computeV(T1, mat->k) / V;
+ Eigen::MatrixXd kmean = HeatTerm::computeMatrix(*e, T1, mat->k) / V;
kmoy11[i] = kmean(0, 0);
kmoy22[i] = kmean(1, 1);
kmoy12[i] = kmean(0, 1);
@@ -378,7 +379,7 @@ void Solver::buildK(Eigen::SparseMatrix<double, Eigen::RowMajor> &K2, Eigen::Map
return;
//std::cout << "processing element #" << e->no << "\n";
- Eigen::MatrixXd Ke = e->buildK(T1, mat->k, false);
+ Eigen::MatrixXd Ke = HeatTerm::build(*e, T1, mat->k, false);
// assembly
tbb::spin_mutex::scoped_lock lock(mutex);
@@ -411,7 +412,7 @@ void Solver::buildK(Eigen::SparseMatrix<double, Eigen::RowMajor> &K2, Eigen::Map
[&](Element *e) {
if (e->type() != ELTYPE::TRI3)
return;
- Eigen::MatrixXd Ke = e->buildK(T1, mat->k, true); // bidon
+ Eigen::MatrixXd Ke = HeatTerm::build(*e, T1, mat->k, true); // bidon
});
}
@@ -436,7 +437,7 @@ void Solver::buildK(Eigen::SparseMatrix<double, Eigen::RowMajor> &K2, Eigen::Map
[&](Element *e) {
if (e->type() != ELTYPE::TRI3)
return;
- Eigen::MatrixXd Ke = e->buildK(T1, mat->k, false);
+ Eigen::MatrixXd Ke = HeatTerm::build(*e, T1, mat->k, false);
// assembly
//tbb::spin_mutex::scoped_lock lock(mutex);
@@ -531,8 +532,7 @@ void Solver::buildK(Eigen::SparseMatrix<double, Eigen::RowMajor> &K2, Eigen::Map
//std::cout << "processing element #" << e->no << "\n";
// ** Ce matrix => K matrix
- Fct0C f(1.0);
- Eigen::VectorXd Ce = e->builds(T1, f);
+ Eigen::VectorXd Ce = HeatTerm::build(*e, T1, Fct0C(1.0));
// assembly
tbb::spin_mutex::scoped_lock lock(mutex);
@@ -603,7 +603,7 @@ void Solver::buildqint(Eigen::VectorXd &qint)
return;
//std::cout << "processing element #" << e->no << "\n";
- Eigen::VectorXd qe = e->computeqint(T1, mat->k);
+ Eigen::VectorXd qe = HeatTerm::build2(*e, T1, mat->k);
// assembly
tbb::spin_mutex::scoped_lock lock(mutex);
@@ -654,7 +654,7 @@ void Solver::builds(Eigen::Map<Eigen::VectorXd> &s)
//std::cout << "processing element #" << e->no << "\n";
// ** se vector => s vector
- Eigen::VectorXd se = e->builds(T1, src->f);
+ Eigen::VectorXd se = HeatTerm::build(*e, T1, src->f);
// assembly
tbb::spin_mutex::scoped_lock lock(mutex);
@@ -694,7 +694,7 @@ void Solver::buildq(Eigen::Map<Eigen::VectorXd> &q)
//std::cout << "processing element #" << e->no << "\n";
// ** qe vector => q vector
- Eigen::VectorXd qe = e->builds(T1, bnd->f);
+ Eigen::VectorXd qe = HeatTerm::build(*e, T1, bnd->f);
// assembly
tbb::spin_mutex::scoped_lock lock(mutex);
diff --git a/mirrors/src/mirrors.h b/mirrors/src/mirrors.h
index e2462ea35ed2eb7094ab68237f4a47b083755d1d..df88ff4147ada4b5010b8c457800924774301d20 100644
--- a/mirrors/src/mirrors.h
+++ b/mirrors/src/mirrors.h
@@ -46,6 +46,7 @@ class TNeumann;
class TSource;
class MSurface;
class ANSYSSolution;
+class ThermoMecaTerm;
}; // namespace mirrors
#endif //MIRRORS_H
diff --git a/mirrors/src/wSolver.cpp b/mirrors/src/wSolver.cpp
index 2038d523c6a5ff925648d764b0887f4dbd18a54f..ef21f944c25d3136485a14c69951b80b2d782b83 100644
--- a/mirrors/src/wSolver.cpp
+++ b/mirrors/src/wSolver.cpp
@@ -24,19 +24,14 @@
#include "wTNeumann.h"
#include "wTSource.h"
#include "wMSurface.h"
+#include "wThermoMecaTerm.h"
#include "wTag.h"
#include "wMshData.h"
#include "wElement.h"
#include "wNode.h"
#include "wANSYSSolution.h"
#include <map>
-#include "wTag.h"
#include <math.h>
-#include "wHex8.h"
-#include "wQuad4.h"
-#include "wTri3.h"
-#include "wCacheHex8.h"
-#include "wGaussHex8.h"
#include "wResults.h"
#include "wMshExport.h"
#include "wLinearSolver.h"
@@ -180,7 +175,7 @@ void Solver::start(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
- Eigen::MatrixXd K_e = e->buildK_mechanical(H);
+ Eigen::MatrixXd K_e = ThermoMecaTerm::buildK(*e, H);
for (size_t ii = 0; ii < e->nodes.size(); ++ii)
{
@@ -215,7 +210,7 @@ void Solver::start(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
- Eigen::MatrixXd L_e = e->buildL_thermal(K_conductivity);
+ Eigen::MatrixXd L_e = ThermoMecaTerm::buildL(*e, K_conductivity);
for (size_t ii = 0; ii < e->nodes.size(); ++ii)
{
@@ -255,7 +250,7 @@ void Solver::start(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
- Eigen::MatrixXd S_e = e->buildS_thermomechanical(D);
+ Eigen::MatrixXd S_e = ThermoMecaTerm::buildS(*e, D);
for (size_t ii = 0; ii < e->nodes.size(); ++ii)
{
@@ -300,7 +295,7 @@ void Solver::start(MshExport *mshWriter)
if (e->type() != ELTYPE::HEX8 && e->type() != ELTYPE::TETRA4)
continue; //if the element is not a tetrahedron or a hexahedron, it is not treated
- Eigen::MatrixXd N_e = e->buildM();
+ Eigen::MatrixXd N_e = ThermoMecaTerm::buildM(*e);
for (size_t ii = 0; ii < e->nodes.size(); ++ii)
{
@@ -343,7 +338,7 @@ void Solver::start(MshExport *mshWriter)
if (e->type() != ELTYPE::HEX8 && e->type() != ELTYPE::TETRA4)
continue; //if the element is not a triangle or a quadrangle, it is not treated
- Eigen::MatrixXd N_e = e->buildS();
+ Eigen::MatrixXd N_e = ThermoMecaTerm::buildM(*e);
for (size_t ii = 0; ii < e->nodes.size(); ++ii)
{
@@ -370,7 +365,7 @@ void Solver::start(MshExport *mshWriter)
if (e->type() != ELTYPE::HEX8 && e->type() != ELTYPE::TETRA4)
continue; //if the element is not a triangle or a quadrangle, it is not treated
- Eigen::MatrixXd N_e = e->buildS();
+ Eigen::MatrixXd N_e = ThermoMecaTerm::buildM(*e);
for (size_t ii = 0; ii < e->nodes.size(); ++ii)
{
@@ -634,7 +629,7 @@ void Solver::start(MshExport *mshWriter)
if (e->type() != ELTYPE::HEX8 && e->type() != ELTYPE::TETRA4)
continue; //if the element is not a triangle or a quadrangle, it is not treated
- Eigen::MatrixXd N_e = e->buildS();
+ Eigen::MatrixXd N_e = ThermoMecaTerm::buildM(*e);
for (size_t ii = 0; ii < e->nodes.size(); ++ii)
{
@@ -748,10 +743,8 @@ void Solver::start(MshExport *mshWriter)
for (auto e : current_tag.elems)
{
- if (e->type() == ELTYPE::HEX8)
- e->strain_stress(Strain, Sigma, H, u);
- else if (e->type() == ELTYPE::TETRA4)
- e->strain_stress(Strain, Sigma, H, u);
+ if (e->type() == ELTYPE::HEX8 || e->type() == ELTYPE::TETRA4)
+ ThermoMecaTerm::strain_stress(*e, H, u, Strain, Sigma);
else
continue;
diff --git a/mirrors/src/wThermoMecaTerm.cpp b/mirrors/src/wThermoMecaTerm.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d8a005e514e73d069390e4e06f0a6c1803aa6447
--- /dev/null
+++ b/mirrors/src/wThermoMecaTerm.cpp
@@ -0,0 +1,179 @@
+/*
+ * 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 k 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 "wThermoMecaTerm.h"
+
+#include "wElement.h"
+#include "wCache.h"
+#include "wGauss.h"
+#include "wMem.h"
+#include "wNode.h"
+
+using namespace mirrors;
+using namespace tbox;
+
+/**
+ * @brief Build volume/surface mass matrix for the thermo-mechanical equation on one element
+ */
+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());
+ for (size_t k = 0; k < gauss.getN(); ++k)
+ {
+ // Shape functions
+ Eigen::VectorXd const &ff = cache.getSf(k);
+ M += ff * ff.transpose() * gauss.getW(k) * mem.getDetJ(k);
+ }
+ return M;
+}
+
+/**
+ * @brief Build mechanical matrix for the thermo-mechanical equation on one element
+ */
+Eigen::MatrixXd ThermoMecaTerm::buildK(tbox::Element const &e, Eigen::MatrixXd const &H)
+{
+ // 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
+ Eigen::MatrixXd K = Eigen::MatrixXd::Zero(3 * ns, 3 * ns); // 3 dimensions
+ for (size_t k = 0; k < gauss.getN(); ++k)
+ {
+ // Jacobian inverse
+ Eigen::Matrix3d const &invJ = mem.getJinv(k);
+ // Fill B matrix (shape functions)
+ Eigen::MatrixXd const &dff = cache.getDsf(k);
+ Eigen::MatrixXd B = Eigen::MatrixXd::Zero(6, 3 * ns);
+ for (size_t i = 0; i < ns; ++i)
+ {
+ // Compute [Ni,x Ni,y, Ni_z]
+ Eigen::Vector3d dN = invJ * dff.col(i);
+ 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.getW(k) * mem.getDetJ(k);
+ }
+ return K;
+}
+
+/**
+ * @brief Build thermal matrix for the thermo-mechanical equation on one element
+ */
+Eigen::MatrixXd ThermoMecaTerm::buildL(Element const &e, Eigen::MatrixXd const &C)
+{
+ // 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);
+ Eigen::MatrixXd const &dff = cache.getDsf(k);
+
+ // [AC] I refactored this from an entry by entry assembly to a full matrix multiplication,
+ // but as nothing tests this method, I could not check the results.
+ L += (C * J * dff).transpose() * J * dff * gauss.getW(k) * detJ;
+ }
+ return L;
+}
+
+/**
+ * @brief Build thermo-mechanical matrix for the thermo-mechanical equation on one element
+ */
+Eigen::MatrixXd ThermoMecaTerm::buildS(Element const &e, Eigen::MatrixXd const &D)
+{
+ // 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
+ Eigen::MatrixXd S = Eigen::MatrixXd::Zero(3 * ns, ns); // 3 dimensions
+ 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);
+ Eigen::MatrixXd const &dff = cache.getDsf(k);
+
+ // S
+ for (size_t i = 0; i < ns; ++i)
+ S.block(i * 3, 0, 3, ns) += D * J * dff.col(i) * cache.getSf(k).transpose() * gauss.getW(k) * detJ;
+ }
+ return S;
+}
+
+/**
+ * @brief Compute strain and stresses on one element
+ */
+void ThermoMecaTerm::strain_stress(Element const &e, Eigen::MatrixXd const &H, std::vector<double> const &u, Eigen::MatrixXd &Strain, Eigen::MatrixXd &Stress)
+{
+ // 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
+ Eigen::MatrixXd u_e(3, ns);
+ // store the solution
+ for (size_t i = 0; i < ns; ++i)
+ {
+ u_e(0, i) = u[3 * (e.nodes[i]->row)];
+ u_e(1, i) = u[3 * (e.nodes[i]->row) + 1];
+ u_e(2, i) = u[3 * (e.nodes[i]->row) + 2];
+ }
+ for (size_t k = 0; k < gauss.getN(); ++k)
+ {
+ // Jacobian and strain
+ Eigen::Matrix3d const &J = mem.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)
+ du_d_Phi += dff.col(i) * u_e.col(i).transpose();
+ Strain += ((J * du_d_Phi) + (J * du_d_Phi).transpose()) / ns / 2;
+ }
+ // stress
+ for (size_t i = 0; i < 3; ++i)
+ for (size_t j = 0; j < 3; ++j)
+ for (size_t k = 0; k < 3; ++k)
+ for (size_t l = 0; l < 3; ++l)
+ Stress(i, j) += H(i * 3 + j, k * 3 + l) * Strain(k, l);
+}
diff --git a/mirrors/src/wThermoMecaTerm.h b/mirrors/src/wThermoMecaTerm.h
new file mode 100644
index 0000000000000000000000000000000000000000..c02674b5be674cdb0958a8e9c114b57875f8a6ce
--- /dev/null
+++ b/mirrors/src/wThermoMecaTerm.h
@@ -0,0 +1,45 @@
+/*
+ * 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 WTHERMOMECATERM_H
+#define WTHERMOMECATERM_H
+
+#include "mirrors.h"
+
+#include <vector>
+#include <Eigen/Dense>
+
+namespace mirrors
+{
+
+/**
+ * @brief Formulation of thermo-mechanical terms
+ */
+class MIRRORS_API ThermoMecaTerm
+{
+public:
+ // Mass matrix
+ static Eigen::MatrixXd buildM(tbox::Element const &e);
+ // Thermo-mechanical matrices
+ static Eigen::MatrixXd buildK(tbox::Element const &e, Eigen::MatrixXd const &H);
+ static Eigen::MatrixXd buildL(tbox::Element const &e, Eigen::MatrixXd const &C);
+ static Eigen::MatrixXd buildS(tbox::Element const &e, Eigen::MatrixXd const &D);
+ // Strain/stress computation
+ static void strain_stress(tbox::Element const &e, Eigen::MatrixXd const &H, std::vector<double> const &u, Eigen::MatrixXd &Strain, Eigen::MatrixXd &Stress);
+};
+
+} // namespace mirrors
+#endif //WTHERMOMECATERM_H
diff --git a/tbox/src/wElement.cpp b/tbox/src/wElement.cpp
index 303f6efad44a566e5b885c2d1c6da4d4627f5e52..aa8fba025e73710a0ced04b034a157f469686070 100644
--- a/tbox/src/wElement.cpp
+++ b/tbox/src/wElement.cpp
@@ -57,31 +57,8 @@ operator<<(std::ostream &out, ELTYPE const &obj)
return out;
}
-TBOX_API std::ostream &
-operator<<(std::ostream &out, MATTYPE const &obj)
-{
- switch (obj)
- {
- case MATTYPE::MAT_K:
- out << "MAT_K";
- break;
- case MATTYPE::MAT_M:
- out << "MAT_M";
- break;
- case MATTYPE::MAT_S:
- out << "MAT_S";
- break;
- default:
- out << "UNKNOWN-" << static_cast<int>(obj);
- break;
- }
- return out;
-}
-
} // namespace tbox
-// ------------------------------------------------------------------------------------
-
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),
@@ -96,85 +73,9 @@ Element::Element(int n, Tag *_ptag, Tag *_etag,
p->elems.push_back(this);
}
-void Element::write(std::ostream &out) const
-{
- out << type() << " #" << no << ":";
-
- out << " nodes=[ ";
- for (auto n : nodes)
- out << "#" << n->no << ' ';
- out << "]";
-
- if (ptag)
- out << " ptag #" << ptag->no;
- if (etag)
- out << " etag #" << etag->no;
- if (!parts.empty())
- {
- out << " parts=[ ";
- for (auto p : parts)
- out << "#" << p->no << ' ';
- out << "]";
- }
-}
-
-/**
- * @brief Test fct - calculate K and print it
- */
-void Element::buildK_test()
-{
- getVMem().update(true);
- std::vector<double> u;
- Eigen::MatrixXd K = buildK(u, Fct0C(1.0));
- std::cout << "K=" << K << '\n';
-}
-
-/**
- * @brief Test fct - calculate M and print it
- */
-void Element::buildM_test()
-{
- getVMem().update(true);
- Eigen::MatrixXd M = buildM();
- std::cout << "M=" << M << '\n';
-}
-
-/**
- * @brief Test fct - calculate S and print it
- */
-void Element::buildS_test()
-{
- getVMem().update(true);
- Eigen::MatrixXd S = buildS();
- std::cout << "S=" << S << '\n';
-}
-
-/**
- * @brief Test fct - calculate s and print it
- */
-void Element::builds_test()
-{
- getVMem().update(true);
- std::vector<double> u;
- Eigen::VectorXd s = builds(u, Fct0C(1.0));
- std::cout << "s=" << s << '\n';
-}
-
-/**
- * @brief Return the values extrapolated at all the element nodes
- * based on the values at Gauss points
- * @authors Kim Liegeois
- */
-Eigen::MatrixXd Element::gp2node(Eigen::MatrixXd &values_at_gp)
+Mem &Element::getVMem() const
{
- Cache &cache = getVCache();
- size_t ngp = cache.getVGauss().getN();
-
- Eigen::MatrixXd A(ngp, nodes.size());
- for (size_t i = 0; i < ngp; ++i)
- A.row(i) = cache.getSf(i).transpose();
-
- return (A.transpose() * A).inverse() * A.transpose() * values_at_gp; // nodes.size() * 18
+ throw std::runtime_error("Element::getVMem not implemented\n");
}
double Element::buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const
@@ -193,37 +94,23 @@ void Element::buildGradientJ(size_t k, std::vector<double> &dDetJ) const
{
throw std::runtime_error("Element::buildGradientJ not implemented\n");
}
-std::vector<double> Element::buildGradientV(int dim) const
-{
- throw std::runtime_error("Element::buildGradientV not implemented\n");
-}
-Eigen::MatrixXd Element::buildK(std::vector<double> const &u, Fct0 const &fct)
-{
- throw std::runtime_error("Element::buildK not implemented\n");
-}
-Eigen::MatrixXd Element::buildK(std::vector<double> const &u, Fct2 const &fct, bool fake)
-{
- throw std::runtime_error("Element::buildK not implemented\n");
-}
-//---------------------
-Eigen::MatrixXd Element::buildK_mechanical(Eigen::MatrixXd const &H)
+
+double Element::computeV() const
{
- throw std::runtime_error("Element::buildK_mechanical not implemented\n");
+ throw std::runtime_error("Element::computeV not implemented\n");
}
-
-Eigen::MatrixXd Element::buildS_thermomechanical(Eigen::MatrixXd const &D)
+std::vector<double> Element::buildGradientV(int dim) const
{
- throw std::runtime_error("Element::buildS_thermomechanical not implemented\n");
+ throw std::runtime_error("Element::buildGradientV not implemented\n");
}
-Eigen::MatrixXd Element::buildL_thermal(Eigen::MatrixXd const &C)
+std::vector<Eigen::Vector3d> Element::computeTangents() const
{
- throw std::runtime_error("Element::buildL_thermial not implemented\n");
+ throw std::runtime_error("Element::computeTangents not implemented\n");
}
-
-void Element::strain_stress(Eigen::MatrixXd &Strain, Eigen::MatrixXd &Stress, Eigen::MatrixXd const &H, std::vector<double> const &u)
+std::vector<Eigen::Vector3d> Element::buildTangents(size_t k) const
{
- throw std::runtime_error("Element::strain_stress not implemented\n");
+ throw std::runtime_error("Element::buildTangents not implemented\n");
}
Eigen::Vector3d Element::normal() const
@@ -235,67 +122,45 @@ std::vector<Eigen::Vector3d> Element::gradientNormal() const
throw std::runtime_error("Element::gradientNormal not implemented\n");
}
-std::vector<Eigen::Vector3d> Element::computeTangents() const
-{
- throw std::runtime_error("Element::computeTangents not implemented\n");
-}
-std::vector<Eigen::Vector3d> Element::buildTangents(size_t k) const
+/**
+ * @brief Return the values extrapolated at all the element nodes
+ * based on the values at Gauss points
+ */
+Eigen::MatrixXd Element::gp2node(Eigen::MatrixXd &values_at_gp) const
{
- throw std::runtime_error("Element::buildTangents not implemented\n");
-}
+ Cache &cache = getVCache();
+ size_t ngp = cache.getVGauss().getN();
-//---------------------
-Eigen::MatrixXd Element::buildM()
-{
- throw std::runtime_error("Element::buildM not implemented\n");
-}
-Eigen::MatrixXd Element::buildS()
-{
- throw std::runtime_error("Element::buildS not implemented\n");
-}
-Eigen::VectorXd Element::builds(std::vector<double> const &u, Fct0 const &f)
-{
- throw std::runtime_error("Element::builds not implemented\n");
-}
-Eigen::VectorXd Element::builds(std::vector<double> const &u, Fct1 const &f)
-{
- throw std::runtime_error("Element::builds not implemented\n");
-}
-Eigen::MatrixXd Element::build(MATTYPE type)
-{
- throw std::runtime_error("Element::build not implemented\n");
-}
-Eigen::VectorXd Element::computeGrad(std::vector<double> const &u, size_t k) const
-{
- throw std::runtime_error("Element::computeGrad not implemented\n");
-}
-double Element::computeGrad2(std::vector<double> const &u, size_t k) const
-{
- throw std::runtime_error("Element::computeGrad2 not implemented\n");
-}
-Eigen::VectorXd Element::computeqV(std::vector<double> const &u, Fct2 const &fct)
-{
- throw std::runtime_error("Element::computeqV not implemented\n");
-}
-Eigen::VectorXd Element::computeqint(std::vector<double> const &u, Fct2 const &fct)
-{
- throw std::runtime_error("Element::computeqint not implemented\n");
-}
-double Element::computeV(std::vector<double> const &u, Fct0 const &fct)
-{
- throw std::runtime_error("Element::computeV not implemented\n");
-}
-Eigen::MatrixXd Element::computeV(std::vector<double> const &u, Fct2 const &fct)
-{
- throw std::runtime_error("Element::computeV not implemented\n");
+ Eigen::MatrixXd A(ngp, nodes.size());
+ for (size_t i = 0; i < ngp; ++i)
+ A.row(i) = cache.getSf(i).transpose();
+
+ return (A.transpose() * A).inverse() * A.transpose() * values_at_gp; // nodes.size() * 18
}
-double Element::computeV()
+Eigen::VectorXd Element::computeGradient(std::vector<double> const &u, size_t k) const
{
- return computeV(std::vector<double>(), Fct0C(1.0));
+ throw std::runtime_error("Element::computeGradient not implemented\n");
}
-Mem &Element::getVMem() const
+void Element::write(std::ostream &out) const
{
- throw std::runtime_error("Element::getVMem not implemented\n");
+ out << type() << " #" << no << ":";
+
+ out << " nodes=[ ";
+ for (auto n : nodes)
+ out << "#" << n->no << ' ';
+ out << "]";
+
+ if (ptag)
+ out << " ptag #" << ptag->no;
+ if (etag)
+ out << " etag #" << etag->no;
+ if (!parts.empty())
+ {
+ out << " parts=[ ";
+ for (auto p : parts)
+ out << "#" << p->no << ' ';
+ out << "]";
+ }
}
diff --git a/tbox/src/wElement.h b/tbox/src/wElement.h
index 068f1bece6d8864928828750312b1d058d4bd37d..72c4dd9b92680ff801c2583ab7c8fa8ef41eeccc 100644
--- a/tbox/src/wElement.h
+++ b/tbox/src/wElement.h
@@ -43,16 +43,6 @@ enum class ELTYPE
TBOX_API std::ostream &operator<<(std::ostream &out, ELTYPE const &obj);
#endif
-enum class MATTYPE
-{
- MAT_K,
- MAT_M,
- MAT_S
-};
-#ifndef SWIG
-TBOX_API std::ostream &operator<<(std::ostream &out, MATTYPE const &obj);
-#endif
-
/**
* @brief Base finite element class
* @authors Romain Boman, Kim Liegeois, Adrien Crovato
@@ -77,48 +67,30 @@ public:
return ELTYPE::UNKNOWN;
}
- // test functions
- void buildK_test();
- void buildM_test();
- void buildS_test();
- void builds_test();
-
#ifndef SWIG
- virtual Eigen::MatrixXd gp2node(Eigen::MatrixXd &values_at_gp);
+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;
- virtual std::vector<double> buildGradientV(int dim) const;
- //--------------------
- virtual void write(std::ostream &out) const override;
- virtual Eigen::MatrixXd buildK(std::vector<double> const &u, Fct0 const &fct);
- virtual Eigen::MatrixXd buildK(std::vector<double> const &u, Fct2 const &fct, bool fake);
- //---------------------
- virtual Eigen::MatrixXd buildK_mechanical(Eigen::MatrixXd const &H);
- virtual Eigen::MatrixXd buildS_thermomechanical(Eigen::MatrixXd const &D);
- virtual Eigen::MatrixXd buildL_thermal(Eigen::MatrixXd const &C);
- virtual void strain_stress(Eigen::MatrixXd &Strain, Eigen::MatrixXd &Stress, Eigen::MatrixXd const &H, std::vector<double> const &u);
- virtual Eigen::Vector3d normal() const;
- virtual std::vector<Eigen::Vector3d> gradientNormal() const;
- virtual std::vector<Eigen::Vector3d> computeTangents() const;
- virtual std::vector<Eigen::Vector3d> buildTangents(size_t k) const;
- //---------------------
- virtual Eigen::MatrixXd buildM();
- virtual Eigen::MatrixXd buildS();
- virtual Eigen::VectorXd builds(std::vector<double> const &u, Fct0 const &f);
- virtual Eigen::VectorXd builds(std::vector<double> const &u, Fct1 const &f);
- virtual Eigen::MatrixXd build(MATTYPE type);
- virtual Eigen::VectorXd computeGrad(std::vector<double> const &u, size_t k) const;
- virtual double computeGrad2(std::vector<double> const &u, size_t k) const;
- virtual Eigen::VectorXd computeqV(std::vector<double> const &u, Fct2 const &fct);
- virtual Eigen::VectorXd computeqint(std::vector<double> const &u, Fct2 const &fct);
- virtual double computeV(std::vector<double> const &u, Fct0 const &fct);
- virtual Eigen::MatrixXd computeV(std::vector<double> const &u, Fct2 const &fct);
- virtual double computeV();
+public:
+ // 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;
+ // Utilities
+ virtual Eigen::VectorXd computeGradient(std::vector<double> const &u, size_t k) const;
+ virtual Eigen::MatrixXd gp2node(Eigen::MatrixXd &values_at_gp) const;
+ virtual void write(std::ostream &out) const override;
#endif
};
diff --git a/tbox/src/wFct2.cpp b/tbox/src/wFct2.cpp
index 4b83977ab0b2b536087ca74f0f56d235a52e08d1..2bf6ac7adb8454b9a666494f23b02a375b73fa31 100644
--- a/tbox/src/wFct2.cpp
+++ b/tbox/src/wFct2.cpp
@@ -106,7 +106,7 @@ void Fct2UdU::eval(Element const &e,
upg += ff(i) * u[e.nodes[i]->row];
//gradU au pg => gradupg
- Eigen::VectorXd gradupg = e.computeGrad(u, k);
+ Eigen::VectorXd gradupg = e.computeGradient(u, k);
this->eval(e, k, upg, gradupg, out, fake);
}
diff --git a/tbox/src/wHex8.cpp b/tbox/src/wHex8.cpp
index da1332f572e05d3ca27238dc4604c04b581ce29a..ae4c523af2c6a7ad95f3fa13b3667711ec455eac 100644
--- a/tbox/src/wHex8.cpp
+++ b/tbox/src/wHex8.cpp
@@ -20,7 +20,6 @@
#include "wCacheHex8.h"
#include "wMemHex8.h"
#include "wNode.h"
-#include "wFct0.h"
using namespace tbox;
@@ -70,179 +69,9 @@ double Hex8::buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const
}
/**
- * @brief Compute stiffness matrix for scalar field
- *
- * K_ij(8,8) = int f*dN_j*dN_i dV
- */
-Eigen::MatrixXd Hex8::buildK(std::vector<double> const &u, Fct0 const &fct)
-{
- CacheHex8 &cache = getCache(order);
- GaussHex8 &gauss = cache.gauss;
- MemHex8 &mem = getMem();
-
- // Gauss integration
- Eigen::MatrixXd K = Eigen::Matrix<double, 8, 8>::Zero();
- for (size_t k = 0; k < gauss.getN(); ++k)
- {
- // Jacobian inverse and shape functions
- Eigen::Matrix3d const &J = mem.getJinv(k);
- Eigen::MatrixXd const &dff = cache.getDsf(k);
-
- // Elementary stiffness matrix
- K += (J * dff).transpose() * J * dff * fct.eval(*this, u, k) * gauss.getW(k) * mem.getDetJ(k);
- }
- return K;
-}
-
-/**
- * @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)
- */
-Eigen::MatrixXd Hex8::buildK_mechanical(Eigen::MatrixXd const &H)
-{
- CacheHex8 &cache = getCache(order);
- GaussHex8 &gauss = cache.gauss;
- MemHex8 &mem = getMem();
-
- Eigen::MatrixXd K = Eigen::MatrixXd::Zero(24, 24);
- for (size_t k = 0; k < gauss.getN(); ++k)
- {
- // Jacobian inverse
- Eigen::Matrix3d const &invJ = mem.getJinv(k);
- // Fill B matrix (shape functions)
- Eigen::MatrixXd const &dff = cache.getDsf(k);
- 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]
- Eigen::Vector3d dN = invJ * dff.col(i);
- 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.getW(k) * mem.getDetJ(k);
- }
- return K;
-}
-
-/**
- * @brief Compute thermomecanical matrix
- */
-Eigen::MatrixXd Hex8::buildS_thermomechanical(Eigen::MatrixXd const &D)
-{
- CacheHex8 &cache = getCache(order);
- GaussHex8 &gauss = cache.gauss;
- MemHex8 &mem = getMem();
-
- Eigen::MatrixXd S_e = Eigen::MatrixXd::Zero(24, 8);
- for (size_t a = 0; a < gauss.getN(); ++a)
- {
- // Jacobian and shape functions
- double detJ = mem.getDetJ(a);
- Eigen::Matrix3d const &J = mem.getJinv(a);
- Eigen::MatrixXd const &dff = cache.getDsf(a);
-
- // S
- for (size_t i = 0; i < nodes.size(); ++i)
- S_e.block(i * 3, 0, 3, nodes.size()) += D * J * dff.col(i) * cache.getSf(a).transpose() * gauss.getW(a) * detJ;
- }
- return S_e;
-}
-
-/**
- * @brief Compute thermal matrix
+ * @brief Compute the element volume
*/
-Eigen::MatrixXd Hex8::buildL_thermal(Eigen::MatrixXd const &C)
-{
- CacheHex8 &cache = getCache(order);
- GaussHex8 &gauss = cache.gauss;
- MemHex8 &mem = getMem();
-
- // Gauss integration
- Eigen::MatrixXd L_e = Eigen::Matrix<double, 8, 8>::Zero();
- for (size_t a = 0; a < gauss.getN(); ++a)
- {
- // Jacobian
- double detJ = mem.getDetJ(a);
- Eigen::Matrix3d const &J = mem.getJinv(a);
- Eigen::MatrixXd const &dff = cache.getDsf(a);
-
- // [AC] I refactored this from an entry by entry assembly to a full matrix multiplication,
- // but as nothing tests this method, I could not check the results.
- L_e += (C * J * dff).transpose() * J * dff * gauss.getW(a) * detJ;
- }
- return L_e;
-}
-
-void Hex8::strain_stress(Eigen::MatrixXd &Strain, Eigen::MatrixXd &Stress, Eigen::MatrixXd const &H, std::vector<double> const &u)
-{
- CacheHex8 &cache = getCache(order);
- GaussHex8 &gauss = cache.gauss;
- MemHex8 &mem = getMem();
-
- Strain = Stress = Eigen::Matrix3d::Zero();
- Eigen::MatrixXd u_e(3, nodes.size());
- for (size_t i = 0; i < nodes.size(); ++i)
- {
- u_e(0, i) = u[3 * (nodes[i]->row)];
- u_e(1, i) = u[3 * (nodes[i]->row) + 1];
- u_e(2, i) = u[3 * (nodes[i]->row) + 2];
- }
- for (size_t a = 0; a < gauss.getN(); ++a)
- {
- // Jacobian and strain
- Eigen::Matrix3d const &J = mem.getJinv(a);
- Eigen::MatrixXd const &dff = cache.getDsf(a);
- Eigen::Matrix3d du_d_Phi = Eigen::Matrix3d::Zero();
- for (size_t i = 0; i < nodes.size(); ++i)
- du_d_Phi += dff.col(i) * u_e.col(i).transpose();
- Strain += ((J * du_d_Phi) + (J * du_d_Phi).transpose()) / nodes.size() / 2;
- }
- // stress
- for (size_t i = 0; i < 3; ++i)
- for (size_t j = 0; j < 3; ++j)
- for (size_t k = 0; k < 3; ++k)
- for (size_t l = 0; l < 3; ++l)
- 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
- */
-Eigen::MatrixXd Hex8::buildM()
-{
- CacheHex8 &cache = getCache(order);
- GaussHex8 &gauss = cache.gauss;
- MemHex8 &mem = getMem();
-
- // Gauss integration
- Eigen::MatrixXd M = Eigen::Matrix<double, 8, 8>::Zero();
- for (size_t k = 0; k < gauss.getN(); ++k)
- {
- // Shape functions
- Eigen::VectorXd const &ff = cache.getSf(k);
- M += ff * ff.transpose() * gauss.getW(k) * mem.getDetJ(k);
- }
- return M;
-}
-
-/**
- * @brief Compute intergal of scalar function
- *
- * I = int f dV
- */
-double Hex8::computeV(std::vector<double> const &u, Fct0 const &fct)
+double Hex8::computeV() const
{
CacheHex8 &cache = getCache(order);
GaussHex8 &gauss = cache.gauss;
@@ -251,14 +80,6 @@ double Hex8::computeV(std::vector<double> const &u, Fct0 const &fct)
// Gauss integration
double V = 0.0;
for (size_t k = 0; k < gauss.getN(); ++k)
- V += fct.eval(*this, u, k) * gauss.getW(k) * mem.getDetJ(k);
+ V += gauss.getW(k) * mem.getDetJ(k);
return V;
}
-
-/**
- * @brief Return the element volume
- */
-double Hex8::computeV()
-{
- return getMem().getVol();
-}
diff --git a/tbox/src/wHex8.h b/tbox/src/wHex8.h
index 1c563df033a1ca4381e156159434978d06f87461..bf39aaa5123a926d59a7182cb9d8da7b329b0a32 100644
--- a/tbox/src/wHex8.h
+++ b/tbox/src/wHex8.h
@@ -33,21 +33,17 @@ 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;
- virtual double buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const override;
- virtual Eigen::MatrixXd buildK(std::vector<double> const &u, Fct0 const &fct) override;
- virtual Eigen::MatrixXd buildM() override;
- virtual Eigen::MatrixXd buildK_mechanical(Eigen::MatrixXd const &H) override;
- virtual Eigen::MatrixXd buildS_thermomechanical(Eigen::MatrixXd const &D) override;
- virtual Eigen::MatrixXd buildL_thermal(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;
+protected:
+ virtual double buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const override;
- virtual double computeV(std::vector<double> const &u, Fct0 const &fct) override;
- virtual double computeV() override;
+public: // @todo
+ virtual double computeV() const override;
#endif
public:
diff --git a/tbox/src/wLine2.cpp b/tbox/src/wLine2.cpp
index 9e34ca2267eda0fc1777103aa84a6d80b1f8db93..85dce33b026da3146b4f99bda70d84fd5bd29a1d 100644
--- a/tbox/src/wLine2.cpp
+++ b/tbox/src/wLine2.cpp
@@ -20,8 +20,6 @@
#include "wCacheLine2.h"
#include "wMemLine2.h"
#include "wNode.h"
-#include "wFct0.h"
-#include "wFct1.h"
using namespace tbox;
@@ -52,59 +50,6 @@ Cache &Line2::getVCache() const
return getCache(order);
}
-/**
- * @brief Compute element tangent vectors
- */
-std::vector<Eigen::Vector3d> Line2::computeTangents() const
-{
- return std::vector<Eigen::Vector3d>{nodes[1]->pos - nodes[0]->pos};
-}
-
-/**
- * @brief Build element tangent vectors using shape functions at Gauss point k
- */
-std::vector<Eigen::Vector3d> Line2::buildTangents(size_t k) const
-{
- Eigen::MatrixXd const &dff = getCache(order).getDsf(k);
-
- std::vector<Eigen::Vector3d> t(1, Eigen::Vector3d::Zero());
- size_t i = 0;
- for (auto it = nodes.begin(); it != nodes.end(); ++it, ++i)
- {
- t[0] += (*it)->pos * dff(0, i);
- }
- return t;
-}
-
-/**
- * @brief Compute element unit normal vector
- */
-Eigen::Vector3d Line2::normal() const
-{
- Eigen::Vector3d t = this->computeTangents()[0].normalized();
- return 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
-{
- // Compute normal and norm
- Eigen::Vector3d t = this->computeTangents()[0];
- Eigen::Vector3d n = Eigen::Vector3d(-t(1), t(0), 0.);
- double nn = n.norm();
-
- // 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;
-}
-
/**
* @brief Compute Jacobian determinant at Gauss point k
*
@@ -141,7 +86,22 @@ void Line2::buildGradientJ(size_t k, std::vector<double> &dDetJ) const
}
/**
- * @brief Compute volume gradient with respect to each node coordinate
+ * @brief Compute the element length
+ */
+double Line2::computeV() const
+{
+ CacheLine2 &cache = getCache(order);
+ GaussLine2 &gauss = cache.gauss;
+ MemLine2 &mem = getMem();
+
+ // Gauss integration
+ double V = 0.0;
+ for (size_t k = 0; k < gauss.getN(); ++k)
+ V += gauss.getW(k) * mem.getDetJ(k);
+ return V;
+}
+/**
+ * @brief Compute length gradient with respect to each node coordinate
*
* dV[i] = w_k * dDetJ[i]_k
*/
@@ -163,74 +123,54 @@ std::vector<double> Line2::buildGradientV(int dim) const
}
/**
- * @brief Compute flux vector
- *
- * s_i(2) = int f * N_i dV
+ * @brief Compute element tangent vectors
*/
-Eigen::VectorXd Line2::builds(std::vector<double> const &u, Fct0 const &f)
+std::vector<Eigen::Vector3d> Line2::computeTangents() const
{
- CacheLine2 &cache = getCache(order);
- GaussLine2 &gauss = cache.gauss;
- MemLine2 &mem = getMem();
-
- // Gauss integration
- Eigen::VectorXd s = Eigen::Vector2d::Zero();
- for (size_t k = 0; k < gauss.getN(); ++k)
- s += cache.getSf(k) * f.eval(*this, u, k) * gauss.getW(k) * mem.getDetJ(k);
- return s;
+ return std::vector<Eigen::Vector3d>{nodes[1]->pos - nodes[0]->pos};
}
/**
- * @brief Compute flux vector
- *
- * s_i(2) = int f_j*n_j * N_i dV
+ * @brief Build element tangent vectors using shape functions at Gauss point k
*/
-Eigen::VectorXd Line2::builds(std::vector<double> const &u, Fct1 const &f)
+std::vector<Eigen::Vector3d> Line2::buildTangents(size_t k) const
{
- CacheLine2 &cache = getCache(order);
- GaussLine2 &gauss = cache.gauss;
- MemLine2 &mem = getMem();
+ Eigen::MatrixXd const &dff = getCache(order).getDsf(k);
- // Gauss integration
- Eigen::VectorXd s = Eigen::Vector2d::Zero();
- for (size_t k = 0; k < gauss.getN(); ++k)
- s += cache.getSf(k) * (f.eval(*this, u, k).dot(normal())) * gauss.getW(k) * mem.getDetJ(k);
- return s;
+ std::vector<Eigen::Vector3d> t(1, Eigen::Vector3d::Zero());
+ size_t i = 0;
+ for (auto it = nodes.begin(); it != nodes.end(); ++it, ++i)
+ {
+ t[0] += (*it)->pos * dff(0, i);
+ }
+ return t;
}
/**
- * @brief Compute square of gradient along element
- *
- * grad = (Delta_U / Delta_L)^2
+ * @brief Compute element unit normal vector
*/
-double Line2::computeGrad2(std::vector<double> const &u, size_t k) const
+Eigen::Vector3d Line2::normal() const
{
- double dUdL = (u[nodes[1]->row] - u[nodes[0]->row]) / getMem().getVol(); // V = length of the element
- return dUdL * dUdL;
+ Eigen::Vector3d t = this->computeTangents()[0].normalized();
+ return Eigen::Vector3d(-t(1), t(0), 0.);
}
/**
- * @brief Compute intergal of scalar function
- *
- * I = int f dV
+ * @brief Compute gradient of element unit normal vector with respect to each node coordinate
*/
-double Line2::computeV(std::vector<double> const &u, Fct0 const &fct)
+std::vector<Eigen::Vector3d> Line2::gradientNormal() const
{
- CacheLine2 &cache = getCache(order);
- GaussLine2 &gauss = cache.gauss;
- MemLine2 &mem = getMem();
-
- // Gauss integration
- double V = 0.0;
- for (size_t k = 0; k < gauss.getN(); ++k)
- V += fct.eval(*this, u, k) * gauss.getW(k) * mem.getDetJ(k);
- return V;
-}
+ // Compute normal and norm
+ Eigen::Vector3d t = this->computeTangents()[0];
+ Eigen::Vector3d n = Eigen::Vector3d(-t(1), t(0), 0.);
+ double nn = n.norm();
-/**
- * @brief Return the element length
- */
-double Line2::computeV()
-{
- return getMem().getVol();
+ // 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;
}
diff --git a/tbox/src/wLine2.h b/tbox/src/wLine2.h
index 4e033f3e3e902153a0012ea4098c4a47f1a2a9cf..9e6e5f0de9dc05c7a676de610640515aec8b1baa 100644
--- a/tbox/src/wLine2.h
+++ b/tbox/src/wLine2.h
@@ -33,20 +33,21 @@ 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;
- virtual std::vector<double> buildGradientV(int dim) const override;
-
- virtual Eigen::VectorXd builds(std::vector<double> const &u, Fct0 const &f) override;
- virtual Eigen::VectorXd builds(std::vector<double> const &u, Fct1 const &f) override;
-
- virtual double computeGrad2(std::vector<double> const &u, size_t k) const override;
- virtual double computeV(std::vector<double> const &u, Fct0 const &fct) override;
- virtual double computeV() 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 std::vector<Eigen::Vector3d> computeTangents() const override;
virtual std::vector<Eigen::Vector3d> buildTangents(size_t k) const override;
#endif
@@ -59,8 +60,6 @@ public:
#ifndef SWIG
virtual Mem &getVMem() const override;
virtual Cache &getVCache() const override;
- virtual Eigen::Vector3d normal() const override;
- virtual std::vector<Eigen::Vector3d> gradientNormal() const override;
#endif
};
diff --git a/tbox/src/wMemHex8.cpp b/tbox/src/wMemHex8.cpp
index aa8398b478c62edfa7b3d58dd50e9a0e0f444744..c815f2025241d25ffaf9c35f40d079999f1c0fdd 100644
--- a/tbox/src/wMemHex8.cpp
+++ b/tbox/src/wMemHex8.cpp
@@ -44,7 +44,7 @@ void MemHex8::update(bool sameDim)
cg /= e.nodes.size();
// Memory is now valid and the volume can be computed
valid = true;
- vol = e.Element::computeV();
+ vol = e.computeV();
// Update gradient if they are needed
if (validG)
diff --git a/tbox/src/wMemLine2.cpp b/tbox/src/wMemLine2.cpp
index 074a6bc6e5b619427bf3978ba19de9a9fc0557b4..8e0c8022ec509f26c9e55bb14802cbb62368c011 100644
--- a/tbox/src/wMemLine2.cpp
+++ b/tbox/src/wMemLine2.cpp
@@ -41,7 +41,7 @@ void MemLine2::update(bool sameDim)
cg /= e.nodes.size();
// Memory is now valid and the volume can be computed
valid = true;
- vol = e.Element::computeV();
+ vol = e.computeV();
// Update gradient if they are needed
if (validG)
diff --git a/tbox/src/wMemQuad4.cpp b/tbox/src/wMemQuad4.cpp
index ad205b1844451cbf14a53c739464d207fbcc8240..50ee3b51eeb5e3fceed1e16ac27ccf2f6e4b2b63 100644
--- a/tbox/src/wMemQuad4.cpp
+++ b/tbox/src/wMemQuad4.cpp
@@ -51,7 +51,7 @@ void MemQuad4::update(bool sameDim)
cg /= e.nodes.size();
// Memory is now valid and the volume can be computed
valid = true;
- vol = e.Element::computeV();
+ vol = e.computeV();
// Update gradient if they are needed
if (validG)
diff --git a/tbox/src/wMemTetra4.cpp b/tbox/src/wMemTetra4.cpp
index 349ca6498ffd7a64d7f5ba5904f317242b8c6991..5be98b0c5ab2aeccd9dd2b28ad58899c27588df2 100644
--- a/tbox/src/wMemTetra4.cpp
+++ b/tbox/src/wMemTetra4.cpp
@@ -44,7 +44,7 @@ void MemTetra4::update(bool sameDim)
cg /= e.nodes.size();
// Memory is now valid and the volume can be computed
valid = true;
- vol = e.Element::computeV();
+ vol = e.computeV();
// Update gradient if they are needed
if (validG)
diff --git a/tbox/src/wMemTri3.cpp b/tbox/src/wMemTri3.cpp
index 4f10a2b806c26c428ead252a1b7ae4eee49aca05..707e7d55838171ad620a1c20f85f5f1b7547aaaf 100644
--- a/tbox/src/wMemTri3.cpp
+++ b/tbox/src/wMemTri3.cpp
@@ -51,7 +51,7 @@ void MemTri3::update(bool sameDim)
cg /= e.nodes.size();
// Memory is now valid and the volume can be computed
valid = true;
- vol = e.Element::computeV();
+ vol = e.computeV();
// Update gradient if they are needed
if (validG)
diff --git a/tbox/src/wMshDeform.cpp b/tbox/src/wMshDeform.cpp
index dd87d8eddd6c32ef9d5e5aab8fe74faba7ed5d86..6eacfebaa80367dea2cc0be5c5976513f3141a54 100644
--- a/tbox/src/wMshDeform.cpp
+++ b/tbox/src/wMshDeform.cpp
@@ -21,6 +21,8 @@
#include "wTag.h"
#include "wNode.h"
#include "wElement.h"
+#include "wCache.h"
+#include "wGauss.h"
#include "wMem.h"
#include "wGmres.h"
@@ -227,6 +229,68 @@ void MshDeform::savePos()
initPos[nDim * i + m] = movBndNodes[i]->pos(m);
}
+/**
+ * @brief Compute mechanical stiffness matrix for one element
+ *
+ * K_ij(6,6) = int eps_k H_ijkl eps_l dV = B(6,3)H(3,3)B(3,6)
+ * K_ij(12,12) = int eps_k H_ijkl eps_l dV = B(12,6)H(6,6)B(6,12)
+ * @todo should be in a dedicated static class
+ */
+Eigen::MatrixXd MshDeform::buildK(tbox::Element const &e, Eigen::MatrixXd const &H)
+{
+ // 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
+ Eigen::MatrixXd K = Eigen::MatrixXd::Zero(nDim * ns, nDim * ns);
+ for (size_t k = 0; k < gauss.getN(); ++k)
+ {
+ // Jacobian inverse
+ Eigen::MatrixXd const &iJ = mem.getJinv(k);
+ // Fill B matrix (shape functions)
+ Eigen::MatrixXd const &dff = cache.getDsf(k);
+ Eigen::MatrixXd B;
+ if (nDim == 2)
+ {
+ B = Eigen::MatrixXd::Zero(3, 2 * ns);
+ for (size_t i = 0; i < ns; ++i)
+ {
+ // Compute [Ni,x Ni,y]
+ Eigen::VectorXd dN = iJ * dff.col(i);
+ 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);
+ }
+ }
+ else
+ {
+ B = Eigen::MatrixXd::Zero(6, 3 * ns);
+ for (size_t i = 0; i < ns; ++i)
+ {
+ // Compute [Ni,x Ni,y Ni,z]
+ Eigen::VectorXd dN = iJ * dff.col(i);
+ 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.getW(k) * mem.getDetJ(k);
+ }
+ return K;
+}
+
/**
* @brief Build the mesh deformation matrix using linear elasticity law
*
@@ -280,7 +344,7 @@ void MshDeform::build(Eigen::SparseMatrix<double, Eigen::RowMajor> &K)
}
// Elementary stiffness matrix
- Eigen::MatrixXd Ke = e->buildK_mechanical(H);
+ Eigen::MatrixXd Ke = this->buildK(*e, H);
// Assembly
// tbb::spin_mutex::scoped_lock lock(mutex);
diff --git a/tbox/src/wMshDeform.h b/tbox/src/wMshDeform.h
index a4a30036d1a649c1165eeed91413c01a7a9f21e3..4547383ae817b8c326fcbcbdbd7301c2ec2e2f76 100644
--- a/tbox/src/wMshDeform.h
+++ b/tbox/src/wMshDeform.h
@@ -58,6 +58,7 @@ private:
std::vector<double> initPos; ///< initial positions of nodes before mesh deformation
void uniqueSort(std::vector<Node *> &nodes);
+ Eigen::MatrixXd buildK(tbox::Element const &e, Eigen::MatrixXd const &H);
public:
std::shared_ptr<MshData> msh;
diff --git a/tbox/src/wQuad4.cpp b/tbox/src/wQuad4.cpp
index c44ab416db48fcde6b60c8115a704aa32c25b04e..d3981e8b5325e0dfd4b6b50447a65b2de5ae17fa 100644
--- a/tbox/src/wQuad4.cpp
+++ b/tbox/src/wQuad4.cpp
@@ -20,7 +20,6 @@
#include "wCacheQuad4.h"
#include "wMemQuad4.h"
#include "wNode.h"
-#include "wFct0.h"
using namespace tbox;
@@ -51,44 +50,6 @@ Cache &Quad4::getVCache() const
return getCache(order);
}
-/**
- * @brief Compute element tangent vectors
- */
-std::vector<Eigen::Vector3d> Quad4::computeTangents() const
-{
- Eigen::Vector3d x1 = nodes[0]->pos;
- Eigen::Vector3d x3 = nodes[2]->pos;
- Eigen::Vector3d x4 = nodes[3]->pos;
-
- return std::vector<Eigen::Vector3d>{(x3 - x1), (x4 - x1)};
-}
-
-/**
- * @brief Build element tangent vectors using shape functions at Gauss point k
- */
-std::vector<Eigen::Vector3d> Quad4::buildTangents(size_t k) const
-{
- Eigen::MatrixXd const &dff = getCache(order).getDsf(k);
-
- std::vector<Eigen::Vector3d> t(2, Eigen::Vector3d::Zero());
- size_t i = 0;
- for (auto it = nodes.begin(); it != nodes.end(); ++it, ++i)
- {
- t[0] += (*it)->pos * dff(0, i);
- t[1] += (*it)->pos * dff(1, i);
- }
- return t;
-}
-
-/**
- * @brief Compute element unit normal vector
- */
-Eigen::Vector3d Quad4::normal() const
-{
- std::vector<Eigen::Vector3d> ts = this->computeTangents();
- return ts[0].cross(ts[1]).normalized();
-}
-
/**
* @brief Compute Jacobian determinant at Gauss point k
*
@@ -118,68 +79,55 @@ double Quad4::buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const
}
/**
- * @brief Compute mass matrix
- *
- * S_ij(4,4) = int N_i*N_j dV
+ * @brief Return the element surface
*/
-Eigen::MatrixXd Quad4::buildS()
+double Quad4::computeV() const
{
CacheQuad4 &cache = getCache(order);
GaussQuad4 &gauss = cache.gauss;
MemQuad4 &mem = getMem();
// Gauss integration
- Eigen::MatrixXd S = Eigen::Matrix4d::Zero();
+ double V = 0.0;
for (size_t k = 0; k < gauss.getN(); ++k)
- {
- // Shape functions
- Eigen::VectorXd const &ff = cache.getSf(k);
-
- // Mass matrix
- S += ff * ff.transpose() * gauss.getW(k) * mem.getDetJ(k);
- }
- return S;
+ V += gauss.getW(k) * mem.getDetJ(k);
+ return V;
}
/**
- * @brief API for build methods
+ * @brief Compute element tangent vectors
*/
-Eigen::MatrixXd Quad4::build(MATTYPE type)
+std::vector<Eigen::Vector3d> Quad4::computeTangents() const
{
- switch (type)
- {
- case MATTYPE::MAT_S:
- return buildS();
- break;
- default:
- std::stringstream out;
- out << "Quad4 cannot build " << type;
- throw std::runtime_error(out.str());
- }
+ Eigen::Vector3d x1 = nodes[0]->pos;
+ Eigen::Vector3d x3 = nodes[2]->pos;
+ Eigen::Vector3d x4 = nodes[3]->pos;
+
+ return std::vector<Eigen::Vector3d>{(x3 - x1), (x4 - x1)};
}
/**
- * @brief Compute intergal of scalar function
- *
- * I = int f dV
+ * @brief Build element tangent vectors using shape functions at Gauss point k
*/
-double Quad4::computeV(std::vector<double> const &u, Fct0 const &fct)
+std::vector<Eigen::Vector3d> Quad4::buildTangents(size_t k) const
{
- CacheQuad4 &cache = getCache(order);
- GaussQuad4 &gauss = cache.gauss;
- MemQuad4 &mem = getMem();
+ Eigen::MatrixXd const &dff = getCache(order).getDsf(k);
- // Gauss integration
- double V = 0.0;
- for (size_t k = 0; k < gauss.getN(); ++k)
- V += fct.eval(*this, u, k) * gauss.getW(k) * mem.getDetJ(k);
- return V;
+ std::vector<Eigen::Vector3d> t(2, Eigen::Vector3d::Zero());
+ size_t i = 0;
+ for (auto it = nodes.begin(); it != nodes.end(); ++it, ++i)
+ {
+ t[0] += (*it)->pos * dff(0, i);
+ t[1] += (*it)->pos * dff(1, i);
+ }
+ return t;
}
/**
- * @brief Return the element surface
+ * @brief Compute element unit normal vector
*/
-double Quad4::computeV()
+Eigen::Vector3d Quad4::normal() const
{
- return getMem().getVol();
-}
\ No newline at end of file
+ std::vector<Eigen::Vector3d> ts = this->computeTangents();
+ return ts[0].cross(ts[1]).normalized();
+}
diff --git a/tbox/src/wQuad4.h b/tbox/src/wQuad4.h
index fa8d1060fc08d652b1df520a6d46ce4e9c37f2f0..0e3975b50c0d2323a7bb2f1c2bde585bc747c375 100644
--- a/tbox/src/wQuad4.h
+++ b/tbox/src/wQuad4.h
@@ -33,20 +33,21 @@ 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;
- virtual Eigen::MatrixXd build(MATTYPE type) override;
- virtual Eigen::MatrixXd buildS() override;
-
- virtual double computeV(std::vector<double> const &u, Fct0 const &fct) override;
- virtual double computeV() override;
-
+public: // @todo
+ virtual double computeV() const 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;
#endif
public:
@@ -57,7 +58,6 @@ public:
#ifndef SWIG
virtual Cache &getVCache() const override;
virtual Mem &getVMem() const override;
- virtual Eigen::Vector3d normal() const override;
#endif
};
diff --git a/tbox/src/wTetra4.cpp b/tbox/src/wTetra4.cpp
index e488be5814bee278b29ed023bd3b1f6afbd2e811..e7b8e5e3716cb6209013558bcc63c23cbd2b460a 100644
--- a/tbox/src/wTetra4.cpp
+++ b/tbox/src/wTetra4.cpp
@@ -20,7 +20,6 @@
#include "wCacheTetra4.h"
#include "wMemTetra4.h"
#include "wNode.h"
-#include "wFct0.h"
using namespace tbox;
@@ -96,204 +95,48 @@ void Tetra4::buildGradientJ(size_t k, std::vector<double> &dDetJ, std::vector<Ei
for (size_t i = 0; i < dDetJ.size(); ++i)
dDetJ[i] = detJ * (iJ * dJ[i]).trace();
}
-/**
- * @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
-{
- GaussTetra4 &gauss = getCache(order).gauss;
- MemTetra4 &mem = getMem();
-
- // Gradient of volume
- std::vector<double> dVol(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
- for (size_t i = 0; i < dDetJ.size(); ++i)
- dVol[i] += w * dDetJ[i];
- }
- return dVol;
-}
/**
- * @brief Compute mass matrix
- *
- * M_ij(4,4) = int N_i*N_j dV
- */
-Eigen::MatrixXd Tetra4::buildM()
-{
- CacheTetra4 &cache = getCache(order);
- GaussTetra4 &gauss = cache.gauss;
- MemTetra4 &mem = getMem();
-
- // Gauss integration
- Eigen::MatrixXd M = Eigen::Matrix4d::Zero();
- for (size_t k = 0; k < gauss.getN(); ++k)
- {
- // Shape functions
- Eigen::VectorXd const &ff = cache.getSf(k);
- M += ff * ff.transpose() * gauss.getW(k) * mem.getDetJ(k);
- }
- return M;
-}
-
-/**
- * @brief Compute stiffness matrix for scalar field
- *
- * K_ij(4,4) = int f*dN_j*dN_i dV
+ * @brief Return the element volume
*/
-Eigen::MatrixXd Tetra4::buildK(std::vector<double> const &u, Fct0 const &fct)
+double Tetra4::computeV() const
{
CacheTetra4 &cache = getCache(order);
GaussTetra4 &gauss = cache.gauss;
MemTetra4 &mem = getMem();
// Gauss integration
- Eigen::MatrixXd K = Eigen::Matrix4d::Zero();
+ double V = 0.0;
for (size_t k = 0; k < gauss.getN(); ++k)
- {
- // Jacobian inverse and shape functions
- Eigen::Matrix3d const &J = mem.getJinv(k);
- Eigen::MatrixXd const &dff = cache.getDsf(k);
-
- // Elementary stiffness matrix
- K += (J * dff).transpose() * J * dff * fct.eval(*this, u, k) * gauss.getW(k) * mem.getDetJ(k);
- }
- return K;
+ V += gauss.getW(k) * mem.getDetJ(k);
+ return V;
}
-
/**
- * @brief Compute linear elasticity stiffness matrix
+ * @brief Compute volume gradient with respect to each node coordinate
*
- * K_ij(12,12) = int eps_k H_ijkl eps_l dV = B(12,6)H(6,6)B(6,12)
+ * dV[i] = w_k * dDetJ[i]_k
*/
-Eigen::MatrixXd Tetra4::buildK_mechanical(Eigen::MatrixXd const &H)
+std::vector<double> Tetra4::buildGradientV(int dim) const
{
- CacheTetra4 &cache = getCache(order);
- GaussTetra4 &gauss = cache.gauss;
+ GaussTetra4 &gauss = getCache(order).gauss;
MemTetra4 &mem = getMem();
- Eigen::MatrixXd K = Eigen::MatrixXd::Zero(12, 12);
+ // Gradient of volume
+ std::vector<double> dVol(4 * dim, 0.); // 4 nodes
for (size_t k = 0; k < gauss.getN(); ++k)
{
- // Jacobian inverse
- Eigen::Matrix3d const &invJ = mem.getJinv(k);
- // Fill B matrix (shape functions)
- Eigen::MatrixXd const &dff = cache.getDsf(k);
- 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]
- Eigen::Vector3d dN = invJ * dff.col(i);
- 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.getW(k) * mem.getDetJ(k);
- }
- return K;
-}
-
-/**
- * @brief Compute thermomechanical mass matrix
- */
-Eigen::MatrixXd Tetra4::buildS_thermomechanical(Eigen::MatrixXd const &D)
-{
- CacheTetra4 &cache = getCache(order);
- GaussTetra4 &gauss = cache.gauss;
- MemTetra4 &mem = getMem();
-
- // Gauss integration
- Eigen::MatrixXd S_e = Eigen::MatrixXd::Zero(12, 4);
- for (size_t a = 0; a < gauss.getN(); ++a)
- {
- // Jacobian and shape functions
- double detJ = mem.getDetJ(a);
- Eigen::Matrix3d const &J = mem.getJinv(a);
- Eigen::MatrixXd const &dff = cache.getDsf(a);
-
- // S
- for (size_t i = 0; i < nodes.size(); ++i)
- S_e.block(i * 3, 0, 3, nodes.size()) += D * J * dff.col(i) * cache.getSf(a).transpose() * gauss.getW(a) * detJ;
- }
- return S_e;
-}
-
-/**
- * @brief Compute thermal matrix
- */
-Eigen::MatrixXd Tetra4::buildL_thermal(Eigen::MatrixXd const &C)
-{
- CacheTetra4 &cache = getCache(order);
- GaussTetra4 &gauss = cache.gauss;
- MemTetra4 &mem = getMem();
-
- // Gauss integration
- Eigen::MatrixXd L_e = Eigen::Matrix4d::Zero();
- for (size_t a = 0; a < gauss.getN(); ++a)
- {
- // Jacobian
- double detJ = mem.getDetJ(a);
- Eigen::Matrix3d const &J = mem.getJinv(a);
- Eigen::MatrixXd const &dff = cache.getDsf(a);
-
- // [AC] I refactored this from an entry by entry assembly to a full matrix multiplication,
- // but as nothing tests this method, I could not check the results.
- L_e += (C * J * dff).transpose() * J * dff * gauss.getW(a) * detJ;
- }
- return L_e;
-}
-
-/**
- * @brief Compute stress and strain
- */
-void Tetra4::strain_stress(Eigen::MatrixXd &Strain, Eigen::MatrixXd &Stress, Eigen::MatrixXd const &H, std::vector<double> const &u)
-{
- CacheTetra4 &cache = getCache(order);
- GaussTetra4 &gauss = cache.gauss;
- MemTetra4 &mem = getMem();
-
- Strain = Stress = Eigen::Matrix3d::Zero();
- Eigen::MatrixXd u_e(3, nodes.size());
- for (size_t i = 0; i < nodes.size(); ++i)
- {
- u_e(0, i) = u[3 * (nodes[i]->row)];
- u_e(1, i) = u[3 * (nodes[i]->row) + 1];
- u_e(2, i) = u[3 * (nodes[i]->row) + 2];
- }
- for (size_t a = 0; a < gauss.getN(); ++a)
- {
- // Jacobian and strain
- Eigen::Matrix3d const &J = mem.getJinv(a);
- Eigen::MatrixXd const &dff = cache.getDsf(a);
- Eigen::Matrix3d du_d_Phi = Eigen::Matrix3d::Zero();
- for (size_t i = 0; i < nodes.size(); ++i)
- du_d_Phi += dff.col(i) * u_e.col(i).transpose();
- Strain += ((J * du_d_Phi) + (J * du_d_Phi).transpose()) / nodes.size() / 2;
+ double w = gauss.getW(k); // GP weight
+ std::vector<double> dDetJ = mem.getGradDetJ(k); // gradient of Jacobian determinant
+ for (size_t i = 0; i < dDetJ.size(); ++i)
+ dVol[i] += w * dDetJ[i];
}
- // stress
- for (size_t i = 0; i < 3; ++i)
- for (size_t j = 0; j < 3; ++j)
- for (size_t k = 0; k < 3; ++k)
- for (size_t l = 0; l < 3; ++l)
- Stress(i, j) += H(i * 3 + j, k * 3 + l) * Strain(k, l);
+ return dVol;
}
/**
* @brief Compute gradient of a scalar field at Gauss point k
*/
-Eigen::VectorXd Tetra4::computeGrad(std::vector<double> const &u, size_t k) const
+Eigen::VectorXd Tetra4::computeGradient(std::vector<double> const &u, size_t k) const
{
// Solution vector
Eigen::Vector4d ue;
@@ -302,37 +145,3 @@ Eigen::VectorXd Tetra4::computeGrad(std::vector<double> const &u, size_t k) cons
// gradient in global axis: inv(J)_ij * d_jN_i * ue_i
return getMem().getJinv(k) * getCache(order).getDsf(k) * ue;
}
-
-/**
- * @brief Compute square of gradient of a scalar field at Gauss point k
- */
-double Tetra4::computeGrad2(std::vector<double> const &u, size_t k) const
-{
- return computeGrad(u, k).squaredNorm();
-}
-
-/**
- * @brief Compute intergal of scalar function
- *
- * I = int f dV
- */
-double Tetra4::computeV(std::vector<double> const &u, Fct0 const &fct)
-{
- CacheTetra4 &cache = getCache(order);
- GaussTetra4 &gauss = cache.gauss;
- MemTetra4 &mem = getMem();
-
- // Gauss integration
- double V = 0.0;
- for (size_t k = 0; k < gauss.getN(); ++k)
- V += fct.eval(*this, u, k) * gauss.getW(k) * mem.getDetJ(k);
- return V;
-}
-
-/**
- * @brief Return the element volume
- */
-double Tetra4::computeV()
-{
- return getMem().getVol();
-}
diff --git a/tbox/src/wTetra4.h b/tbox/src/wTetra4.h
index f524f9c9243994d196274a59768b400e3fb0a6a0..cc2bdb3970b4939322dfad8537dae06162d94798 100644
--- a/tbox/src/wTetra4.h
+++ b/tbox/src/wTetra4.h
@@ -33,25 +33,19 @@ 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;
- virtual std::vector<double> buildGradientV(int dim) const override;
- virtual Eigen::MatrixXd buildM() override;
- virtual Eigen::MatrixXd buildK(std::vector<double> const &u, Fct0 const &fct) override;
- virtual Eigen::MatrixXd buildK_mechanical(Eigen::MatrixXd const &H) override;
-
- virtual Eigen::MatrixXd buildS_thermomechanical(Eigen::MatrixXd const &D) override;
- virtual Eigen::MatrixXd buildL_thermal(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;
-
- virtual Eigen::VectorXd computeGrad(std::vector<double> const &u, size_t k) const override;
- virtual double computeGrad2(std::vector<double> const &u, size_t k) const override;
- virtual double computeV(std::vector<double> const &u, Fct0 const &fct) override;
- virtual double computeV() override;
+public: // @todo
+ virtual double computeV() const override;
+ virtual std::vector<double> buildGradientV(int dim) const override;
#endif
public:
@@ -62,6 +56,7 @@ public:
#ifndef SWIG
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/wTri3.cpp b/tbox/src/wTri3.cpp
index 2be7fa463b9b24f7b95fa9024294bcf767a2a7c4..f145c4ed2bc8b947ed188c6e1c990f7f54e34849 100644
--- a/tbox/src/wTri3.cpp
+++ b/tbox/src/wTri3.cpp
@@ -20,9 +20,6 @@
#include "wCacheTri3.h"
#include "wMemTri3.h"
#include "wNode.h"
-#include "wFct0.h"
-#include "wFct1.h"
-#include "wFct2.h"
using namespace tbox;
@@ -53,69 +50,6 @@ Cache &Tri3::getVCache() const
return getCache(order);
}
-/**
- * @brief Compute element tangent vectors
- */
-std::vector<Eigen::Vector3d> Tri3::computeTangents() const
-{
- Eigen::Vector3d x1 = nodes[0]->pos;
- Eigen::Vector3d x2 = nodes[1]->pos;
- Eigen::Vector3d x3 = nodes[2]->pos;
-
- return std::vector<Eigen::Vector3d>{(x2 - x1), (x3 - x1)};
-}
-
-/**
- * @brief Build element tangent vectors using shape functions at Gauss point k
- */
-std::vector<Eigen::Vector3d> Tri3::buildTangents(size_t k) const
-{
- Eigen::MatrixXd const &dff = getCache(order).getDsf(k);
-
- std::vector<Eigen::Vector3d> t(2, Eigen::Vector3d::Zero());
- size_t i = 0;
- for (auto it = nodes.begin(); it != nodes.end(); ++it, ++i)
- {
- t[0] += (*it)->pos * dff(0, i);
- t[1] += (*it)->pos * dff(1, i);
- }
- return t;
-}
-
-/**
- * @brief Compute element unit normal vector
- */
-Eigen::Vector3d Tri3::normal() const
-{
- std::vector<Eigen::Vector3d> ts = this->computeTangents();
- return 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
-{
- // Compute normal and norm
- std::vector<Eigen::Vector3d> ts = this->computeTangents();
- Eigen::Vector3d n = ts[0].cross(ts[1]);
- double nn = n.norm();
-
- // 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
- 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
- }
- }
- return dn;
-}
-
/**
* @brief Compute Jacobian determinant at Gauss point k
*
@@ -197,6 +131,22 @@ void Tri3::buildGradientJ(size_t k, std::vector<double> &dDetJ, std::vector<Eige
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
+{
+ CacheTri3 &cache = getCache(order);
+ GaussTri3 &gauss = cache.gauss;
+ MemTri3 &mem = getMem();
+
+ // Gauss integration
+ double V = 0.0;
+ for (size_t k = 0; k < gauss.getN(); ++k)
+ V += gauss.getW(k) * mem.getDetJ(k);
+ return V;
+}
/**
* @brief Compute volume gradient with respect to each node coordinate
*
@@ -220,168 +170,72 @@ std::vector<double> Tri3::buildGradientV(int dim) const
}
/**
- * @brief Compute stiffness matrix for scalar field
- *
- * K_ij(3,3) = int f*dN_j*dN_i dV
- */
-Eigen::MatrixXd Tri3::buildK(std::vector<double> const &u, Fct0 const &fct)
-{
- CacheTri3 &cache = getCache(order);
- GaussTri3 &gauss = cache.gauss;
- MemTri3 &mem = getMem();
-
- // Gauss integration
- Eigen::MatrixXd K = Eigen::Matrix3d::Zero();
- for (size_t k = 0; k < gauss.getN(); ++k)
- {
- // Jacobian inverse and shape functions
- Eigen::Matrix2d const &J = mem.getJinv(k);
- Eigen::MatrixXd const &dff = cache.getDsf(k);
-
- // Elementary stiffness matrix
- K += (J * dff).transpose() * J * dff * fct.eval(*this, u, k) * gauss.getW(k) * mem.getDetJ(k);
- }
- return K;
-}
-
-/**
- * @brief Compute stiffness matrix for scalar field
- *
- * K_ij(3,3) = int [f]*dN_j*dN_i dV
+ * @brief Compute element tangent vectors
*/
-Eigen::MatrixXd Tri3::buildK(std::vector<double> const &u, Fct2 const &fct, bool fake)
+std::vector<Eigen::Vector3d> Tri3::computeTangents() const
{
- CacheTri3 &cache = getCache(order);
- GaussTri3 &gauss = cache.gauss;
-
- Eigen::MatrixXd K = Eigen::Matrix3d::Zero();
- if (fake) // fake run - fill MPI job list
- {
- for (size_t k = 0; k < gauss.getN(); ++k)
- {
- Eigen::MatrixXd fk;
- fct.eval(*this, u, k, fk, fake);
- }
- }
- else // true run - use MPI cached results
- {
- // Gauss integration
- MemTri3 &mem = getMem();
- for (size_t k = 0; k < gauss.getN(); ++k)
- {
- // Jacobian inverse, shape functions and function evaluation
- Eigen::Matrix2d const &J = mem.getJinv(k);
- Eigen::MatrixXd const &dff = cache.getDsf(k);
- Eigen::MatrixXd fk;
- fct.eval(*this, u, k, fk, fake);
+ Eigen::Vector3d x1 = nodes[0]->pos;
+ Eigen::Vector3d x2 = nodes[1]->pos;
+ Eigen::Vector3d x3 = nodes[2]->pos;
- // Elementary stiffness matrix
- K += (fk * J * dff).transpose() * J * dff * gauss.getW(k) * mem.getDetJ(k);
- }
- }
- return K;
+ return std::vector<Eigen::Vector3d>{(x2 - x1), (x3 - x1)};
}
/**
- * @brief Compute mechanical stiffness matrix
- *
- * K_ij(6,6) = int eps_k H_ijkl eps_l dV = B(6,3)H(3,3)B(3,6)
+ * @brief Build element tangent vectors using shape functions at Gauss point k
*/
-Eigen::MatrixXd Tri3::buildK_mechanical(Eigen::MatrixXd const &H)
+std::vector<Eigen::Vector3d> Tri3::buildTangents(size_t k) const
{
- CacheTri3 &cache = getCache(order);
- GaussTri3 &gauss = cache.gauss;
- MemTri3 &mem = getMem();
+ Eigen::MatrixXd const &dff = getCache(order).getDsf(k);
- // Gauss integration
- Eigen::MatrixXd K = Eigen::MatrixXd::Zero(6, 6);
- for (size_t k = 0; k < gauss.getN(); ++k)
+ std::vector<Eigen::Vector3d> t(2, Eigen::Vector3d::Zero());
+ size_t i = 0;
+ for (auto it = nodes.begin(); it != nodes.end(); ++it, ++i)
{
- // Jacobian inverse
- Eigen::Matrix2d const &invJ = mem.getJinv(k);
- // Fill B matrix (shape functions)
- Eigen::MatrixXd const &dff = cache.getDsf(k);
- Eigen::MatrixXd B = Eigen::MatrixXd::Zero(3, 6);
- for (size_t i = 0; i < nodes.size(); ++i)
- {
- // Compute [Ni,x Ni,y]
- Eigen::Vector2d dN = invJ * dff.col(i);
- 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 * gauss.getW(k) * mem.getDetJ(k);
+ t[0] += (*it)->pos * dff(0, i);
+ t[1] += (*it)->pos * dff(1, i);
}
- return K;
-}
-
-/**
- * @brief Compute flux vector for a scalar field
- *
- * s_i(3) = int f*N_i dV
- */
-Eigen::VectorXd Tri3::builds(std::vector<double> const &u, Fct0 const &f)
-{
- CacheTri3 &cache = getCache(order);
- GaussTri3 &gauss = cache.gauss;
- MemTri3 &mem = getMem();
-
- // Gauss integration
- Eigen::VectorXd s = Eigen::Vector3d::Zero();
- for (size_t k = 0; k < gauss.getN(); ++k)
- s += cache.getSf(k) * f.eval(*this, u, k) * gauss.getW(k) * mem.getDetJ(k);
- return s;
+ return t;
}
/**
- * @brief Compute flux vector for scalar field
- *
- * s_i(3) = int f_j*n_j * N_i dV
+ * @brief Compute element unit normal vector
*/
-Eigen::VectorXd Tri3::builds(std::vector<double> const &u, Fct1 const &f)
+Eigen::Vector3d Tri3::normal() const
{
- CacheTri3 &cache = getCache(order);
- GaussTri3 &gauss = cache.gauss;
- MemTri3 &mem = getMem();
-
- // Gauss integration
- Eigen::VectorXd s = Eigen::Vector3d::Zero();
- for (size_t k = 0; k < gauss.getN(); ++k)
- s += cache.getSf(k) * (f.eval(*this, u, k).dot(normal())) * gauss.getW(k) * mem.getDetJ(k);
- return s;
+ std::vector<Eigen::Vector3d> ts = this->computeTangents();
+ return ts[0].cross(ts[1]).normalized();
}
/**
- * @brief Compute mass matrix
- *
- * S_ij(3,3) = int N_i*N_j dV
+ * @brief Compute gradient of element unit normal vector with respect to each node coordinate
*/
-Eigen::MatrixXd Tri3::buildS()
+std::vector<Eigen::Vector3d> Tri3::gradientNormal() const
{
- CacheTri3 &cache = getCache(order);
- GaussTri3 &gauss = cache.gauss;
- MemTri3 &mem = getMem();
+ // Compute normal and norm
+ std::vector<Eigen::Vector3d> ts = this->computeTangents();
+ Eigen::Vector3d n = ts[0].cross(ts[1]);
+ double nn = n.norm();
- // Gauss integration
- Eigen::MatrixXd S = Eigen::Matrix3d::Zero();
- for (size_t k = 0; k < gauss.getN(); ++k)
+ // 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
+ for (size_t i = 0; i < 3; ++i)
{
- // Shape functions
- Eigen::VectorXd const &ff = cache.getSf(k);
-
- // Mass matrix
- S += ff * ff.transpose() * gauss.getW(k) * mem.getDetJ(k);
+ 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
+ }
}
- return S;
+ return dn;
}
/**
* @brief Compute gradient of a scalar field at Gauss point k
*/
-Eigen::VectorXd Tri3::computeGrad(std::vector<double> const &u, size_t k) const
+Eigen::VectorXd Tri3::computeGradient(std::vector<double> const &u, size_t k) const
{
// Solution vector
Eigen::Vector3d ue;
@@ -390,112 +244,3 @@ Eigen::VectorXd Tri3::computeGrad(std::vector<double> const &u, size_t k) const
// Gradient in global axis: inv(J)_ij * d_jN_i * ue_i
return getMem().getJinv(k) * getCache(order).getDsf(k) * ue;
}
-
-/**
- * @brief Compute square of gradient of a scalar field at Gauss point k
- */
-double Tri3::computeGrad2(std::vector<double> const &u, size_t k) const
-{
- return computeGrad(u, k).squaredNorm();
-}
-
-/**
- * @brief Compute flux
- *
- * qV_i(3) = int [f] gradu N_i dV
- */
-Eigen::VectorXd Tri3::computeqint(std::vector<double> const &u, Fct2 const &fct)
-{
- CacheTri3 &cache = getCache(order);
- GaussTri3 &gauss = cache.gauss;
- MemTri3 &mem = getMem();
-
- // Gauss integration
- Eigen::VectorXd qV = Eigen::Vector3d::Zero();
- for (size_t k = 0; k < gauss.getN(); ++k)
- {
- // Shape functions, gradient and flux evaluation
- Eigen::MatrixXd const &dff = cache.getDsf(k);
- Eigen::MatrixXd fk(2, 2);
- fct.eval(*this, u, k, fk, false);
-
- // Elementary flux vector
- qV += (fk * computeGrad(u, k)).transpose() * mem.getJinv(k) * dff * gauss.getW(k) * mem.getDetJ(k);
- }
- return qV;
-}
-
-/**
- * @brief Compute flux
- *
- * qV_i(3) = int -[f] gradu dV
- */
-Eigen::VectorXd Tri3::computeqV(std::vector<double> const &u, Fct2 const &fct)
-{
- CacheTri3 &cache = getCache(order);
- GaussTri3 &gauss = cache.gauss;
- MemTri3 &mem = getMem();
-
- // Gauss integration
- Eigen::VectorXd qV = Eigen::Vector2d::Zero();
- for (size_t k = 0; k < gauss.getN(); ++k)
- {
- // Gradient and matrix function evaluation
- Eigen::VectorXd gradk = computeGrad(u, k);
- Eigen::MatrixXd fk(2, 2);
- fct.eval(*this, u, k, fk, false);
-
- qV -= fk * gradk * gauss.getW(k) * mem.getDetJ(k);
- }
- return qV;
-}
-
-/**
- * @brief Compute intergal of scalar function
- *
- * I = int f dV
- */
-double Tri3::computeV(std::vector<double> const &u, Fct0 const &fct)
-{
- CacheTri3 &cache = getCache(order);
- GaussTri3 &gauss = cache.gauss;
- MemTri3 &mem = getMem();
-
- // Gauss integration
- double V = 0.0;
- for (size_t k = 0; k < gauss.getN(); ++k)
- V += fct.eval(*this, u, k) * gauss.getW(k) * mem.getDetJ(k);
- return V;
-}
-
-/**
- * @brief Compute intergal of matrix function
- *
- * I = int [f] dV
- */
-Eigen::MatrixXd Tri3::computeV(std::vector<double> const &u, Fct2 const &fct)
-{
- CacheTri3 &cache = getCache(order);
- GaussTri3 &gauss = cache.gauss;
- MemTri3 &mem = getMem();
-
- // Gauss integration
- Eigen::MatrixXd out = Eigen::Matrix2d::Zero(); /* @todo out should be resized to match fk size */
- for (size_t k = 0; k < gauss.getN(); ++k)
- {
- // Function evaluation
- Eigen::MatrixXd fk;
- fct.eval(*this, u, k, fk, false);
-
- out += fk * gauss.getW(k) * mem.getDetJ(k);
- }
- return out;
-}
-
-/**
- * @brief Return the element surface
- */
-double Tri3::computeV()
-{
- return getMem().getVol();
-}
diff --git a/tbox/src/wTri3.h b/tbox/src/wTri3.h
index bea667503ed240fc0e2ba09a4d214813d8057ccb..bb82eab1bb13448d34c42172d74eb1476c92a0d9 100644
--- a/tbox/src/wTri3.h
+++ b/tbox/src/wTri3.h
@@ -33,30 +33,23 @@ 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;
- virtual std::vector<double> buildGradientV(int dim) const override;
-
- virtual Eigen::MatrixXd buildK(std::vector<double> const &u, Fct0 const &fct) override;
- virtual Eigen::MatrixXd buildK(std::vector<double> const &u, Fct2 const &fct, bool fake) override;
- virtual Eigen::MatrixXd buildK_mechanical(Eigen::MatrixXd const &H) override;
- virtual Eigen::VectorXd builds(std::vector<double> const &u, Fct0 const &f) override;
- virtual Eigen::VectorXd builds(std::vector<double> const &u, Fct1 const &f) override;
- virtual Eigen::MatrixXd buildS() override;
-
- virtual Eigen::VectorXd computeGrad(std::vector<double> const &u, size_t k) const override;
- virtual double computeGrad2(std::vector<double> const &u, size_t k) const override;
- virtual Eigen::VectorXd computeqV(std::vector<double> const &u, Fct2 const &fct) override;
- virtual Eigen::VectorXd computeqint(std::vector<double> const &u, Fct2 const &fct) override;
- virtual double computeV(std::vector<double> const &u, Fct0 const &fct) override;
- virtual Eigen::MatrixXd computeV(std::vector<double> const &u, Fct2 const &fct) override;
- virtual double computeV() 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 std::vector<Eigen::Vector3d> computeTangents() const override;
virtual std::vector<Eigen::Vector3d> buildTangents(size_t k) const override;
#endif
@@ -69,8 +62,7 @@ public:
#ifndef SWIG
virtual Mem &getVMem() const override;
virtual Cache &getVCache() const override;
- virtual Eigen::Vector3d normal() const override;
- virtual std::vector<Eigen::Vector3d> gradientNormal() const override;
+ virtual Eigen::VectorXd computeGradient(std::vector<double> const &u, size_t k) const override;
#endif
};
diff --git a/tbox/tests/fem/elem1s.py b/tbox/tests/fem/elem1s.py
index 6bcca62a71048ed581a859f1597bcb0ed047e9a0..a14f11b9117986a5ce5f4f47db6eccd0eb027019 100644
--- a/tbox/tests/fem/elem1s.py
+++ b/tbox/tests/fem/elem1s.py
@@ -13,10 +13,6 @@ def main(**d):
el = msh.ntags['S1'].elems[0]
print(el)
- #el.buildM_test()
- el.buildK_test()
- el.buildS_test()
- el.builds_test()
# 1 hex
print('='*80+'\nTesting Hex\n'+'='*80)
@@ -24,10 +20,6 @@ def main(**d):
el = msh.ntags['Volume'].elems[0]
print(el)
- el.buildM_test()
- el.buildK_test()
- #el.buildS_test()
- #el.builds_test()
# 1 quad
print('='*80+'\nTesting Quad\n'+'='*80)
@@ -36,10 +28,6 @@ def main(**d):
el = msh.ntags['S1'].elems[0]
print(el)
- #el.buildM_test()
- #el.buildK_test()
- el.buildS_test()
- #el.builds_test()
if __name__ == "__main__":
diff --git a/waves/src/wTimeIntegration.cpp b/waves/src/wTimeIntegration.cpp
index b2603d6755851293d9fe897ff8bf9349667d948a..d958627397747f8a6bd35e424fcfcee28f11d3ed 100644
--- a/waves/src/wTimeIntegration.cpp
+++ b/waves/src/wTimeIntegration.cpp
@@ -24,6 +24,7 @@
#include "wProblem.h"
#include "wMedium.h"
#include "wBoundary.h"
+#include "wWaveTerm.h"
#include "wFct0.h"
#include "wMshExport.h"
#include <typeinfo>
@@ -108,7 +109,7 @@ void TimeIntegration::buildS(Eigen::SparseMatrix<double, Eigen::RowMajor> &S2)
//std::cout << "processing element #" << e->no << "\n";
// ** Se matrix => S vector
- Eigen::MatrixXd Se = e->buildS();
+ Eigen::MatrixXd Se = WaveTerm::buildM(*e);
// assembly
tbb::spin_mutex::scoped_lock lock(mutex);
@@ -165,10 +166,10 @@ void TimeIntegration::buildKM_tbb_lambda(Eigen::SparseMatrix<double, Eigen::RowM
//std::cout << "processing element #" << e->no << "\n";
// ** Me matrix => Md vector
- Eigen::MatrixXd Me = e->buildM();
+ Eigen::MatrixXd Me = WaveTerm::buildM(*e);
// ** Ke matrix => K matrix
- Eigen::MatrixXd Ke = e->buildK(u, Fct0C(1.0));
+ Eigen::MatrixXd Ke = WaveTerm::buildK(*e, u);
// assembly
tbb::spin_mutex::scoped_lock lock(mutex);
@@ -195,7 +196,7 @@ void TimeIntegration::buildKM_tbb_lambda(Eigen::SparseMatrix<double, Eigen::RowM
std::cout << "[cpu] " << chrono1 << '\n';
}
-void TimeIntegration::build(MATTYPE type, Eigen::SparseMatrix<double, Eigen::RowMajor> &A2)
+/*void TimeIntegration::build(MATTYPE type, Eigen::SparseMatrix<double, Eigen::RowMajor> &A2)
{
tbb::spin_mutex mutex;
tbb::task_scheduler_init init(nthreads); // TODO mettre ça ailleurs...
@@ -235,7 +236,7 @@ void TimeIntegration::build(MATTYPE type, Eigen::SparseMatrix<double, Eigen::Row
std::cout << "S (" << A2.rows() << "," << A2.cols() << ") nnz=" << A2.nonZeros() << "\n";
chrono1.read();
std::cout << "[cpu] " << chrono1 << '\n';
-}
+}*/
void TimeIntegration::write(std::ostream &out) const
{
diff --git a/waves/src/wTimeIntegration.h b/waves/src/wTimeIntegration.h
index e34250c2e746a565c45af4d453bd754704e953d4..824bde9b425a29531ce5710d7a23a2713af705a5 100644
--- a/waves/src/wTimeIntegration.h
+++ b/waves/src/wTimeIntegration.h
@@ -78,7 +78,7 @@ public:
void setGUI(DisplayHook &hook) { dhook = &hook; }
- void build(tbox::MATTYPE type, Eigen::SparseMatrix<double, Eigen::RowMajor> &A2);
+ //void build(tbox::MATTYPE type, Eigen::SparseMatrix<double, Eigen::RowMajor> &A2);
void stop() { stopit = true; }
diff --git a/waves/src/wWaveTerm.cpp b/waves/src/wWaveTerm.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e39582ec569b9782b05a1324817b0b2f6de3306e
--- /dev/null
+++ b/waves/src/wWaveTerm.cpp
@@ -0,0 +1,70 @@
+/*
+ * 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 "wWaveTerm.h"
+
+#include "wElement.h"
+#include "wCache.h"
+#include "wGauss.h"
+#include "wMem.h"
+
+using namespace waves;
+using namespace tbox;
+
+/**
+ * @brief Build volume stiffness matrix for the waves equation on one element
+ */
+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::MatrixXd const &dff = cache.getDsf(k);
+
+ // Elementary stiffness matrix
+ K += (J * dff).transpose() * J * dff * gauss.getW(k) * mem.getDetJ(k);
+ }
+ return K;
+}
+
+/**
+ * @brief Build volume/surface mass matrix for the waves equation on one element
+ */
+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());
+ for (size_t k = 0; k < gauss.getN(); ++k)
+ {
+ // Shape functions
+ Eigen::VectorXd const &ff = cache.getSf(k);
+ M += ff * ff.transpose() * gauss.getW(k) * mem.getDetJ(k);
+ }
+ return M;
+}
diff --git a/waves/src/wWaveTerm.h b/waves/src/wWaveTerm.h
new file mode 100644
index 0000000000000000000000000000000000000000..12e60a89f0a9d8c2c3b67e8fa7acb3bba6d17e37
--- /dev/null
+++ b/waves/src/wWaveTerm.h
@@ -0,0 +1,41 @@
+/*
+ * 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 WWAVETERM_H
+#define WWAVETERM_H
+
+#include "waves.h"
+
+#include <vector>
+#include <Eigen/Dense>
+
+namespace waves
+{
+
+/**
+ * @brief Formulation of wave terms
+ */
+class WAVES_API WaveTerm
+{
+public:
+ // Stiffness matrix
+ static Eigen::MatrixXd buildK(tbox::Element const &e, std::vector<double> const &u);
+ // Mass matrix
+ static Eigen::MatrixXd buildM(tbox::Element const &e);
+};
+
+} // namespace waves
+#endif //WWAVETERM_H
diff --git a/waves/src/waves.h b/waves/src/waves.h
index 300ba7cb385088f8fc8d7f133c0fbc271eb0da86..d63fdfd5ba2dc66de0552e04d21e7128d29f8041 100644
--- a/waves/src/waves.h
+++ b/waves/src/waves.h
@@ -40,6 +40,8 @@ class Boundary;
class DisplayHook;
+class WaveTerm;
+
class TimeIntegration;
class ForwardEuler;
class RungeKutta;