feat(thermal): add local shell kernel and time grid

src/fesa/elements/thermal_shell.cpp

@@ -0,0 +1,43 @@
+#include "fesa/elements/thermal_shell.h"
+
+#include <array>
+#include <cmath>
+#include <stdexcept>
+
+namespace fesa {
+namespace {
+constexpr std::array<double, 3> kGauss{{-0.7745966692414834, 0.0, 0.7745966692414834}};
+constexpr std::array<double, 3> kWeight{{5.0 / 9.0, 8.0 / 9.0, 5.0 / 9.0}};
+constexpr std::array<double, 4> kXi{{-1., 1., 1., -1.}};
+constexpr std::array<double, 4> kEta{{-1., -1., 1., 1.}};
+void AddZ(Matrix& z0, Matrix& z1, std::size_t first) {
+  const double a[3][3]={{4,2,-1},{2,16,2},{-1,2,4}};
+  const double b[3][3]={{7,-8,1},{-8,16,-8},{1,-8,7}};
+  for (std::size_t i=0;i<3;++i) for (std::size_t j=0;j<3;++j) { z0(first+i,first+j)+=a[i][j]/30.; z1(first+i,first+j)+=b[i][j]/3.; }
+}
+}  // namespace
+
+ThermalShell::ThermalShell(std::array<std::size_t, 4> node_indices,
+                           std::array<Vector3, 4> coordinates,
+                           double conductivity, double density,
+                           double specific_heat, double thickness)
+    : node_indices_{node_indices} {
+  if (!(conductivity > 0.) || !(density > 0.) || !(specific_heat > 0.) ||
+      !(thickness > 0.)) throw std::invalid_argument("invalid thermal shell material");
+  Matrix m{4U,4U}; Matrix s{4U,4U}; Matrix z0{5U,5U}; Matrix z1{5U,5U};
+  AddZ(z0,z1,0U); AddZ(z0,z1,2U);
+  for (std::size_t x=0;x<3;++x) for (std::size_t y=0;y<3;++y) {
+    const double xi=kGauss[x], eta=kGauss[y]; std::array<double,4> n{}, dx{}, dy{};
+    Vector3 a1{},a2{};
+    for(std::size_t i=0;i<4;++i) { n[i]=.25*(1+kXi[i]*xi)*(1+kEta[i]*eta); dx[i]=.25*kXi[i]*(1+kEta[i]*eta); dy[i]=.25*kEta[i]*(1+kXi[i]*xi); a1=a1+dx[i]*coordinates[i]; a2=a2+dy[i]*coordinates[i]; }
+    const Vector3 w=a1.Cross(a2); const double jac=w.Norm(); const double g11=a1.Dot(a1),g12=a1.Dot(a2),g22=a2.Dot(a2), det=g11*g22-g12*g12;
+    if (!(jac>0.) || !(det>0.) || !std::isfinite(jac)) throw std::invalid_argument("invalid thermal shell geometry");
+    const double wt=kWeight[x]*kWeight[y]*jac;
+    for(std::size_t i=0;i<4;++i) for(std::size_t j=0;j<4;++j) { m(i,j)+=n[i]*n[j]*wt; s(i,j)+=(dx[i]*(g22*dx[j]-g12*dy[j])+dy[i]*(g11*dy[j]-g12*dx[j]))/det*wt; }
+  }
+  for(std::size_t i=0;i<4;++i) for(std::size_t j=0;j<4;++j) for(std::size_t a=0;a<5;++a) for(std::size_t b=0;b<5;++b) { const std::size_t p=5*i+a,q=5*j+b; capacity_(p,q)=density*specific_heat*thickness*.5*m(i,j)*z0(a,b); conductivity_(p,q)=conductivity*(thickness*.5*s(i,j)*z0(a,b)+2./thickness*m(i,j)*z1(a,b)); }
+}
+const std::array<std::size_t,4>& ThermalShell::NodeIndices() const noexcept { return node_indices_; }
+const Matrix& ThermalShell::Conductivity() const noexcept { return conductivity_; }
+const Matrix& ThermalShell::Capacity() const noexcept { return capacity_; }
+}  // namespace fesa