phynexis-deps/mfem
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Pre-update commit containing much temporary code

Browse Source
This commit is contained in:
Stowell, Mark L
2024-05-02 09:40:33 -07:00
parent 1d46fbba57
commit c5be56264c
22 changed files with 3056 additions and 1202 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fem/fe/fe_rt.cpp
+522
View File
@@ -1264,6 +1264,528 @@ void RT_WedgeElement::CalcDivShape(const IntegrationPoint &ip,
}
}
const real_t RT_FuentesPyramidElement::nk[15] =
{ 0,0,-1, 0,-1,0, 1,0,1, 0,1,1, -1,0,0 };
RT_FuentesPyramidElement::RT_FuentesPyramidElement(const int p)
: VectorFiniteElement(3, Geometry::PYRAMID, (p + 1)*(3*p*(p + 2) + 5),
p + 1, H_DIV, FunctionSpace::Pk),
dof2nk(dof)
{
zmax = 0.0;
std::cout << "p " << p << ", order " << order << ", dof " << dof << std::endl;
const real_t *iop = poly1d.OpenPoints(p);
const real_t *icp = poly1d.ClosedPoints(p + 1);
const real_t *bop = poly1d.OpenPoints(p);
std::cout << "have points " << std::endl;
#ifndef MFEM_THREAD_SAFE
tmp1_i.SetSize(p + 2);
tmp1_ij.SetSize(p + 2, p + 2);
tmp2_ij.SetSize(p + 2, dim);
tmp3_ij.SetSize(p + 2, dim);
tmp1_ijk.SetSize(p + 1, p + 1, dim);
tmp2_ijk.SetSize(p + 1, p + 1, dim);
tmp3_ijk.SetSize(p + 1, p + 1, dim);
tmp4_ijk.SetSize(p + 1, p + 2, dim);
tmp5_ijk.SetSize(p + 1, p + 2, dim);
tmp6_ijk.SetSize(p + 2, p + 2, dim);
tmp7_ijk.SetSize(p + 2, p + 2, dim);
u.SetSize(dof, dim);
divu.SetSize(dof);
#else
// Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
Vector tmp1_i(p + 2);
DenseMatrix tmp1_ij(p + 2, p + 2);
DenseMatrix tmp2_ij(p + 2, dim);
DenseMatrix tmp3_ij(p + 2, dim);
DenseTensor tmp1_ijk(p + 1, p + 1, dim);
DenseTensor tmp2_ijk(p + 1, p + 1, dim);
DenseTensor tmp3_ijk(p + 1, p + 1, dim);
DenseTensor tmp4_ijk(p + 1, p + 2, dim);
DenseTensor tmp5_ijk(p + 1, p + 2, dim);
DenseTensor tmp6_ijk(p + 2, p + 2, dim);
DenseTensor tmp7_ijk(p + 2, p + 2, dim);
DenseMatrix u(dof, dim);
#endif
int o = 0;
// quadrilateral face
for (int j = 0; j <= p; j++)
for (int i = 0; i <= p; i++) // (0,1,2,3)
{
Nodes.IntPoint(o).Set3(bop[i], bop[j], 0.);
dof2nk[o++] = 0;
}
// triangular faces
for (int j = 0; j <= p; j++)
for (int i = 0; i + j <= p; i++) // (0,1,4)
{
real_t w = bop[i] + bop[j] + bop[p-i-j];
Nodes.IntPoint(o).Set3(bop[i]/w, 0., bop[j]/w);
dof2nk[o++] = 1;
}
for (int j = 0; j <= p; j++)
for (int i = 0; i + j <= p; i++) // (1,2,4)
{
real_t w = bop[i] + bop[j] + bop[p-i-j];
Nodes.IntPoint(o).Set3(1.-bop[j]/w, bop[i]/w, bop[j]/w);
dof2nk[o++] = 2;
}
for (int j = 0; j <= p; j++)
for (int i = 0; i + j <= p; i++) // (3,4,2)
{
real_t w = bop[i] + bop[j] + bop[p-i-j];
Nodes.IntPoint(o).Set3(bop[j]/w, 1.0-bop[i]/w, bop[i]/w);
dof2nk[o++] = 3;
}
for (int j = 0; j <= p; j++)
for (int i = 0; i + j <= p; i++) // (0,4,3)
{
real_t w = bop[i] + bop[j] + bop[p-i-j];
Nodes.IntPoint(o).Set3(0., bop[j]/w, bop[i]/w);
dof2nk[o++] = 4;
}
// interior
// x-components
for (int k = 0; k <= p; k++)
for (int j = 0; j <= p; j++)
for (int i = 1; i <= p; i++)
{
real_t w = 1.0 - iop[k];
Nodes.IntPoint(o).Set3(icp[i]*w, iop[j]*w, iop[k]);
dof2nk[o++] = 4;
}
// y-components
for (int k = 0; k <= p; k++)
for (int j = 1; j <= p; j++)
for (int i = 0; i <= p; i++)
{
real_t w = 1.0 - iop[k];
Nodes.IntPoint(o).Set3(iop[i]*w, icp[j]*w, iop[k]);
dof2nk[o++] = 1;
}
// z-components
for (int k = 1; k <= p; k++)
for (int j = 0; j <= p; j++)
for (int i = 0; i <= p; i++)
{
real_t w = 1.0 - icp[k];
Nodes.IntPoint(o).Set3(iop[i]*w, iop[j]*w, icp[k]);
dof2nk[o++] = 0;
}
std::cout << "Nodes are set" << std::endl;
DenseMatrix T(dof);
for (int m = 0; m < dof; m++)
{
const IntegrationPoint &ip = Nodes.IntPoint(m);
const Vector nm({nk[3*dof2nk[m]], nk[3*dof2nk[m]+1], nk[3*dof2nk[m]+2]});
std::cout << "calling calcBasis for point " << m << " with normal " << nm(
0) << " " << nm(1) << " " << nm(2) << std::endl;
calcBasis(order, ip, tmp1_i, tmp1_ij, tmp2_ij,
tmp1_ijk, tmp2_ijk, tmp3_ijk, tmp4_ijk, tmp5_ijk, tmp6_ijk,
tmp7_ijk,
tmp3_ij, u);
// std::cout << "u = "; u.Print(std::cout, 10);
u.Mult(nm, T.GetColumn(m));
}
Ti.Factor(T);
}
void RT_FuentesPyramidElement::CalcVShape(const IntegrationPoint &ip,
DenseMatrix &shape) const
{
const int p = order - 1;
#ifdef MFEM_THREAD_SAFE
DenseMatrix u(dof, dim);
#endif
// calcBasis(p, ip, shape_0, shape_1, shape_2, u);
calcBasis(order, ip, tmp1_i, tmp1_ij, tmp2_ij,
tmp1_ijk, tmp2_ijk, tmp3_ijk, tmp4_ijk, tmp5_ijk, tmp6_ijk,
tmp7_ijk, tmp3_ij, u);
Ti.Mult(u, shape);
}
void RT_FuentesPyramidElement::CalcRawVShape(const IntegrationPoint &ip,
DenseMatrix &shape) const
{
const int p = order - 1;
#ifdef MFEM_THREAD_SAFE
// DenseMatrix u(dof, dim);
#endif
// calcBasis(p, ip, shape_0, shape_1, shape_2, u);
calcBasis(order, ip, tmp1_i, tmp1_ij, tmp2_ij,
tmp1_ijk, tmp2_ijk, tmp3_ijk, tmp4_ijk, tmp5_ijk, tmp6_ijk,
tmp7_ijk, tmp3_ij, shape);
// Ti.Mult(u, shape);
}
void RT_FuentesPyramidElement::CalcDivShape(const IntegrationPoint &ip,
Vector &divshape) const
{
const int p = order - 1;
#ifdef MFEM_THREAD_SAFE
// Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
// Vector dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1), dshape_l(p + 1);
Vector divu(dof);
#endif
divu = 0.0;
/*
poly1d.CalcBasis(p, ip.x, shape_x, dshape_x);
poly1d.CalcBasis(p, ip.y, shape_y, dshape_y);
poly1d.CalcBasis(p, ip.z, shape_z, dshape_z);
poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l, dshape_l);
int o = 0;
for (int k = 0; k <= p; k++)
for (int j = 0; j + k <= p; j++)
for (int i = 0; i + j + k <= p; i++)
{
int l = p - i - j - k;
divu(o++) = (dshape_x(i)*shape_l(l) -
shape_x(i)*dshape_l(l))*shape_y(j)*shape_z(k);
divu(o++) = (dshape_y(j)*shape_l(l) -
shape_y(j)*dshape_l(l))*shape_x(i)*shape_z(k);
divu(o++) = (dshape_z(k)*shape_l(l) -
shape_z(k)*dshape_l(l))*shape_x(i)*shape_y(j);
}
for (int j = 0; j <= p; j++)
for (int i = 0; i + j <= p; i++)
{
int k = p - i - j;
divu(o++) =
(shape_x(i) + (ip.x - c)*dshape_x(i))*shape_y(j)*shape_z(k) +
(shape_y(j) + (ip.y - c)*dshape_y(j))*shape_x(i)*shape_z(k) +
(shape_z(k) + (ip.z - c)*dshape_z(k))*shape_x(i)*shape_y(j);
}
*/
Ti.Mult(divu, divshape);
}
void RT_FuentesPyramidElement::calcBasis(const int p,
const IntegrationPoint &ip,
Vector &phi_k,
DenseMatrix &phi_ij,
DenseMatrix &dphi_k,
DenseTensor &VQ_ijk,
DenseTensor &VTa_ijk,
DenseTensor &VTb_ijk,
DenseTensor &E_ijk,
DenseTensor &dE_ijk,
DenseTensor &dphi_ijk,
DenseTensor &VL_ijk,
DenseMatrix &VR_ij,
DenseMatrix &u) const
{
real_t x = ip.x;
real_t y = ip.y;
real_t z = ip.z;
Vector xy({x,y});
real_t mu, muInv;
// Vector muNu(3), dmuNu(3);
if (std::fabs(1.0 - z) < 1e-4)
{
std::cout << "z is close to 1: " << 1.0 - z << std::endl;
}
zmax = std::max(z, zmax);
u = 0.0;
int o = 0;
// Quadrilateral face
if (z < 1.0)
{
V_Q(p, mu01(z, xy, 1), mu01_grad_mu01(z, xy, 1),
mu01(z, xy, 2), mu01_grad_mu01(z, xy, 2),
VQ_ijk);
/*
std::cout << "k = " << 0 << std::endl;
V1_ijk(0).Print(std::cout);
std::cout << "k = " << 1 << std::endl;
V1_ijk(1).Print(std::cout);
std::cout << "k = " << 2 << std::endl;
V1_ijk(2).Print(std::cout);
*/
const real_t muz3 = pow(mu0(z), 3);
for (int j=0; j<p; j++)
for (int i=0; i<p; i++, o++)
for (int k=0; k<3; k++)
{
u(o, k) = muz3 * VQ_ijk(i, j, k);
}
}
// Triangular faces
if (z < 1.0)
{
// (a,b) = (1,2), c = 0
V_T(p, nu012(z, xy, 1), nu012_grad_nu012(z, xy, 1), VTa_ijk);
mu = mu0(z, xy, 2);
if (mu > 0.0)
{
muInv = 1.0 / mu;
//muNu.Set(mu, nu012(z, xy, 1));
//dmuNu.Set(pow(mu, 3), nu012_grad_nu012(z, xy, 1));
V_T(p, lam125(x, y, z), lam125_grad_lam125(x, y, z), VTb_ijk);
}
else
{
//std::cout << "mu02(" << x << "," << y << "," << z << ") <= 0, " << mu << std::endl;
muInv = 1.0;
VTb_ijk = 0.0;
}
for (int j=0; j<p; j++)
for (int i=0; i+j<p; i++, o++)
for (int k=0; k<3; k++)
{
u(o, k) = 0.5 * (mu * VTa_ijk(i, j, k) +
muInv * VTb_ijk(i, j, k));
}
// (a,b) = (1,2), c = 1
mu = mu1(z, xy, 2);
if (mu > 0.0)
{
//std::cout << "mu12(" << x << "," << y << "," << z << ") > 0, " << mu << std::endl;
muInv = 1.0 / mu;
//muNu.Set(mu, nu012(z, xy, 1));
//dmuNu.Set(pow(mu, 3), nu012_grad_nu012(z, xy, 1));
V_T(p, lam435(x, y, z), lam435_grad_lam435(x, y, z), VTb_ijk);
}
else
{
//std::cout << "mu12(" << x << "," << y << "," << z << ") <= 0, " << mu << std::endl;
muInv = 1.0;
VTb_ijk = 0.0;
}
for (int j=0; j<p; j++)
for (int i=0; i+j<p; i++, o++)
for (int k=0; k<3; k++)
{
u(o, k) = 0.5 * (mu * VTa_ijk(i, j, k) +
muInv * VTb_ijk(i, j, k));
}
// (a,b) = (2,1), c = 0
V_T(p, nu012(z, xy, 2), nu012_grad_nu012(z, xy, 2), VTa_ijk);
mu = mu0(z, xy, 1);
if (mu > 0.0)
{
muInv = 1.0 / mu;
//muNu.Set(mu, nu012(z, xy, 2));
//dmuNu.Set(pow(mu, 3), nu012_grad_nu012(z, xy, 2));
//V_T(p, muNu, dmuNu, VTb_ijk);
V_T(p, lam145(x, y, z), lam145_grad_lam145(x, y, z), VTb_ijk);
}
else
{
//std::cout << "mu01(" << x << "," << y << "," << z << ") <= 0, " << mu << std::endl;
muInv = 1.0;
VTb_ijk = 0.0;
}
for (int j=0; j<p; j++)
for (int i=0; i+j<p; i++, o++)
for (int k=0; k<3; k++)
{
u(o, k) = 0.5 * (mu * VTa_ijk(i, j, k) +
muInv * VTb_ijk(i, j, k));
}
// (a,b) = (2,1), c = 1
mu = mu1(z, xy, 1);
muInv = 1.0 / mu;
if (mu > 0.0)
{
//muNu.Set(mu, nu012(z, xy, 2));
//dmuNu.Set(pow(mu, 3), nu012_grad_nu012(z, xy, 2));
//V_T(p, muNu, dmuNu, VTb_ijk);
V_T(p, lam235(x, y, z), lam235_grad_lam235(x, y, z), VTb_ijk);
}
else
{
//std::cout << "mu11(" << x << "," << y << "," << z << ") <= 0, " << mu << std::endl;
muInv = 1.0;
VTb_ijk = 0.0;
}
for (int j=0; j<p; j++)
for (int i=0; i+j<p; i++, o++)
for (int k=0; k<3; k++)
{
u(o, k) = 0.5 * (mu * VTa_ijk(i, j, k) +
muInv * VTb_ijk(i, j, k));
}
}
// Interior
// Family I
if (z < 1.0 && p >= 2)
{
// std::cout << "Calling E_Q at " << x << " " << y << " " << z << std::endl;
E_Q(p, mu01(z, xy, 1), grad_mu01(z, xy, 1),
mu01(z, xy, 2), grad_mu01(z, xy, 2), E_ijk, dE_ijk);
// std::cout << "Calling phi_E at " << x << " " << y << " " << z << std::endl;
phi_E(p, mu01(z), grad_mu01(z), phi_k, dphi_k);
// std::cout << "phi_k(2) " << phi_k(2) << std::endl;
// std::cout << "dphi_k(2, 0) " << dphi_k(2, 0) << std::endl;
// std::cout << "dphi_k(2, 1) " << dphi_k(2, 1) << std::endl;
// std::cout << "dphi_k(2, 2) " << dphi_k(2, 2) << std::endl;
const real_t muz = mu0(z);
const Vector dmuz(grad_mu0(z));
Vector dmuphi(3), E(3), v(3);
for (int k=2; k<=p; k++)
{
dmuphi(0) = muz * dphi_k(k,0) + dmuz(0) * phi_k(k);
dmuphi(1) = muz * dphi_k(k,1) + dmuz(1) * phi_k(k);
dmuphi(2) = muz * dphi_k(k,2) + dmuz(2) * phi_k(k);
for (int j=2; j<=p; j++)
for (int i=0; i<p; i++, o++)
{
E(0) = E_ijk(i,j,0); E(1) = E_ijk(i,j,1); E(2) = E_ijk(i,j,2);
dmuphi.cross3D(E, v);
for (int l=0; l<3; l++)
{
u(o, l) = muz * phi_k(k) * dE_ijk(i,j,l) + v(l);
}
}
}
}
// Family II
if (z < 1.0 && p >= 2)
{
E_Q(p, mu01(z, xy, 2), grad_mu01(z, xy, 2),
mu01(z, xy, 1), grad_mu01(z, xy, 1), E_ijk, dE_ijk);
// phi_E(p, mu01(z), grad_mu01(z), phi_k, dphi_k);
const real_t muz = mu0(z);
const Vector dmuz(grad_mu0(z));
Vector dmuphi(3), E(3), v(3);
for (int k=2; k<=p; k++)
{
dmuphi(0) = muz * dphi_k(k,0) + dmuz(0) * phi_k(k);
dmuphi(1) = muz * dphi_k(k,1) + dmuz(1) * phi_k(k);
dmuphi(2) = muz * dphi_k(k,2) + dmuz(2) * phi_k(k);
for (int j=2; j<=p; j++)
for (int i=0; i<p; i++, o++)
{
E(0) = E_ijk(i,j,0); E(1) = E_ijk(i,j,1); E(2) = E_ijk(i,j,2);
dmuphi.cross3D(E, v);
for (int l=0; l<3; l++)
{
u(o, l) = muz * phi_k(k) * dE_ijk(i,j,l) + v(l);
}
}
}
}
// Family III
if (z < 1.0 && p >= 2)
{
phi_Q(p, mu01(z, xy, 2), grad_mu01(z, xy, 2),
mu01(z, xy, 1), grad_mu01(z, xy, 1), phi_ij, dphi_ijk);
const real_t muz = mu0(z);
const Vector dmuz(grad_mu0(z));
for (int j=2; j<=p; j++)
for (int i=2; i<=p; i++, o++)
{
const int n = std::max(i,j);
const real_t nmu = n * pow(muz, n-1);
u(o, 0) = nmu * (dphi_ijk(i,j,1) * dmuz(2) -
dphi_ijk(i,j,2) * dmuz(1));
u(o, 1) = nmu * (dphi_ijk(i,j,2) * dmuz(0) -
dphi_ijk(i,j,0) * dmuz(2));
u(o, 2) = nmu * (dphi_ijk(i,j,0) * dmuz(1) -
dphi_ijk(i,j,1) * dmuz(0));
}
}
// Family IV
if (z < 1.0 && p >= 2)
{
/*
V_Q(p, mu01(z, xy, 1), mu01_grad_mu01(z, xy, 1),
mu01(z, xy, 2), mu01_grad_mu01(z, xy, 2),
V1_ijk);
*/
phi_E(p, mu01(z), phi_k);
const real_t muz2 = pow(mu0(z), 2);
for (int k=2; k<=p; k++)
for (int j=0; j<p; j++)
for (int i=0; i<p; i++, o++)
for (int l=0; l<3; l++)
{
u(o, l) = muz2 * VQ_ijk(i, j, l) * phi_k(k);
}
}
// Family V
if (z < 1.0 && p >= 2)
{
V_L(p, mu01(z, xy, 1), grad_mu01(z, xy, 1),
mu01(z, xy, 2), grad_mu01(z, xy, 2), mu0(z), grad_mu0(z), VL_ijk);
const real_t muz = mu1(z);
for (int j=2; j<=p; j++)
for (int i=2; i<=p; i++, o++)
{
const int n = std::max(i, j);
const real_t muzi = pow(muz, n-1);
for (int l=0; l<3; l++)
{
u(o, l) = muzi * VL_ijk(i, j, l);
}
}
}
// Family VI
if (z < 1.0 && p >= 2)
{
V_R(p, mu01(z, xy, 1), grad_mu01(z, xy, 1),
mu1(z, xy, 2), grad_mu1(z, xy, 2), mu0(z), grad_mu0(z), VR_ij);
const real_t muz = mu1(z);
for (int i=2; i<=p; i++, o++)
{
const real_t muzi = pow(muz, i-1);
for (int l=0; l<3; l++)
{
u(o, l) = muzi * VR_ij(i, l);
}
}
}
// Family VII
if (z < 1.0 && p >= 2)
{
V_R(p, mu01(z, xy, 2), grad_mu01(z, xy, 2),
mu1(z, xy, 1), grad_mu1(z,xy,1), mu0(z), grad_mu0(z), VR_ij);
const real_t muz = mu1(z);
for (int i=2; i<=p; i++, o++)
{
const real_t muzi = pow(muz, i-1);
for (int l=0; l<3; l++)
{
u(o, l) = muzi * VR_ij(i, l);
}
}
}
}
const real_t RT_R1D_SegmentElement::nk[9] = { 1.,0.,0., 0.,1.,0., 0.,0.,1. };
RT_R1D_SegmentElement::RT_R1D_SegmentElement(const int p,