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

Making lim_dist scalar gridfunc

Browse Source
This commit is contained in:
camierjs
2020-08-26 13:37:28 -07:00
parent e2c2800b43
commit 0ab1717cfc
7 changed files with 61 additions and 40 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fem/kernels.hpp
+22 -1
View File
@@ -117,7 +117,7 @@ MFEM_HOST_DEVICE inline void LoadBGt(const int D1D, const int Q1D,
/// Load 2D input scalar into shared memory
template<int MD1, int NBZ>
MFEM_HOST_DEVICE inline void LoadS(const int e, const int D1D,
MFEM_HOST_DEVICE inline void LoadX(const int e, const int D1D,
const DeviceTensor<3, const double> x,
double sX[NBZ][MD1*MD1])
{
@@ -593,6 +593,27 @@ MFEM_HOST_DEVICE inline void GradXt(const int D1D, const int Q1D,
/// Load 3D scalar input vector into shared memory
template<int MD1>
MFEM_HOST_DEVICE inline void LoadX(const int e, const int D1D,
const DeviceTensor<4, const double> x,
double sm[MD1*MD1*MD1])
{
DeviceCube X(sm, MD1, MD1, MD1);
MFEM_FOREACH_THREAD(dz,z,D1D)
{
MFEM_FOREACH_THREAD(dy,y,D1D)
{
MFEM_FOREACH_THREAD(dx,x,D1D)
{
X(dx,dy,dz) = x(dx,dy,dz,e);
}
}
}
MFEM_SYNC_THREAD;
}
/// Load 3D scalar input vector into shared memory, with comp
template<int MD1>
MFEM_HOST_DEVICE inline void LoadX(const int e, const int D1D, const int c,
const DeviceTensor<5, const double> x,
double sm[MD1*MD1*MD1])
fem/tmop_pa_h2s_c0.cpp
+6 -6
View File
@@ -39,7 +39,7 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_C0_2D,
const auto C0 = const_c0 ?
Reshape(c0_.Read(), 1, 1, 1) :
Reshape(c0_.Read(), Q1D, Q1D, NE);
const auto LD = Reshape(lim_dist.Read(), D1D, D1D, DIM, NE);
const auto LD = Reshape(lim_dist.Read(), D1D, D1D, NE);
const auto J = Reshape(j_.Read(), DIM, DIM, Q1D, Q1D, NE);
const auto W = Reshape(w_.Read(), Q1D, Q1D);
const auto b = Reshape(b_.Read(), Q1D, D1D);
@@ -56,9 +56,9 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_C0_2D,
MFEM_SHARED double B[MQ1*MD1];
MFEM_SHARED double XY[2][NBZ][MD1*MD1];
MFEM_SHARED double DQ[2][NBZ][MD1*MQ1];
MFEM_SHARED double QQ[2][NBZ][MQ1*MQ1];
MFEM_SHARED double XY[NBZ][MD1*MD1];
MFEM_SHARED double DQ[NBZ][MD1*MQ1];
MFEM_SHARED double QQ[NBZ][MQ1*MQ1];
kernels::LoadX<MD1,NBZ>(e,D1D,LD,XY);
kernels::LoadB<MD1,MQ1>(D1D,Q1D,b,B);
@@ -75,9 +75,9 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_C0_2D,
const double coeff0 = const_c0 ? C0(0,0,0) : C0(qx,qy,e);
const double weight_m = weight * lim_normal * coeff0;
double D[2];
double D;
kernels::PullEval<MQ1,NBZ>(qx,qy,QQ,D);
const double dist = D[0]; // GetValues, default comp set to 0
const double dist = D; // GetValues, default comp set to 0
// lim_func->Eval_d2(p1, p0, d_vals(q), grad_grad);
// d2.Diag(1.0 / (dist * dist), x.Size());
fem/tmop_pa_h3s_c0.cpp
+7 -7
View File
@@ -38,7 +38,7 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_Kernel_C0_3D,
const auto C0 = const_c0 ?
Reshape(c0_.Read(), 1, 1, 1, 1) :
Reshape(c0_.Read(), Q1D, Q1D, Q1D, NE);
const auto LD = Reshape(lim_dist.Read(), D1D, D1D, D1D, DIM, NE);
const auto LD = Reshape(lim_dist.Read(), D1D, D1D, D1D, NE);
const auto J = Reshape(j_.Read(), DIM, DIM, Q1D, Q1D, Q1D, NE);
const auto W = Reshape(w_.Read(), Q1D, Q1D, Q1D);
const auto b = Reshape(b_.Read(), Q1D, D1D);
@@ -54,10 +54,10 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_Kernel_C0_3D,
constexpr int MD1 = T_D1D ? T_D1D : T_MAX;
MFEM_SHARED double B[MQ1*MD1];
MFEM_SHARED double DDD[3][MD1*MD1*MD1];
MFEM_SHARED double DDQ[3][MD1*MD1*MQ1];
MFEM_SHARED double DQQ[3][MD1*MQ1*MQ1];
MFEM_SHARED double QQQ[3][MQ1*MQ1*MQ1];
MFEM_SHARED double DDD[MD1*MD1*MD1];
MFEM_SHARED double DDQ[MD1*MD1*MQ1];
MFEM_SHARED double DQQ[MD1*MQ1*MQ1];
MFEM_SHARED double QQQ[MQ1*MQ1*MQ1];
kernels::LoadX<MD1>(e,D1D,LD,DDD);
kernels::LoadB<MD1,MQ1>(D1D,Q1D,b,B);
@@ -77,9 +77,9 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_Kernel_C0_3D,
const double coeff0 = const_c0 ? C0(0,0,0,0) : C0(qx,qy,qz,e);
const double weight_m = weight * lim_normal * coeff0;
double D[3];
double D;
kernels::PullEval<MQ1>(qx,qy,qz,QQQ,D);
const double dist = D[0]; // GetValues, default comp set to 0
const double dist = D; // GetValues, default comp set to 0
// lim_func->Eval_d2(p1, p0, d_vals(q), grad_grad);
// d2.Diag(1.0 / (dist * dist), x.Size());
fem/tmop_pa_p2_c0.cpp
+6 -6
View File
@@ -43,7 +43,7 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_2D,
const auto C0 = const_c0 ?
Reshape(c0_.Read(), 1, 1, 1) :
Reshape(c0_.Read(), Q1D, Q1D, NE);
const auto LD = Reshape(lim_dist.Read(), D1D, D1D, DIM, NE);
const auto LD = Reshape(lim_dist.Read(), D1D, D1D, NE);
const auto J = Reshape(j_.Read(), DIM, DIM, Q1D, Q1D, NE);
const auto b = Reshape(b_.Read(), Q1D, D1D);
const auto W = Reshape(w_.Read(), Q1D, Q1D);
@@ -62,9 +62,9 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_2D,
MFEM_SHARED double B[MQ1*MD1];
MFEM_SHARED double XY[2][NBZ][MD1*MD1];
MFEM_SHARED double DQ[2][NBZ][MD1*MQ1];
MFEM_SHARED double QQ[2][NBZ][MQ1*MQ1];
MFEM_SHARED double XY[NBZ][MD1*MD1];
MFEM_SHARED double DQ[NBZ][MD1*MQ1];
MFEM_SHARED double QQ[NBZ][MQ1*MQ1];
MFEM_SHARED double XY0[2][NBZ][MD1*MD1];
MFEM_SHARED double DQ0[2][NBZ][MD1*MQ1];
@@ -97,13 +97,13 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_2D,
const double detJtr = kernels::Det<2>(Jtr);
const double weight = W(qx,qy) * detJtr;
double ld[2], p0[2], p1[2];
double ld, p0[2], p1[2];
const double coeff0 = const_c0 ? C0(0,0,0) : C0(qx,qy,e);
kernels::PullEval<MQ1,NBZ>(qx,qy,QQ,ld);
kernels::PullEval<MQ1,NBZ>(qx,qy,QQ0,p0);
kernels::PullEval<MQ1,NBZ>(qx,qy,QQ1,p1);
const double dist = ld[0]; // GetValues, default comp set to 0
const double dist = ld; // GetValues, default comp set to 0
double d1[2];
// Eval_d1
fem/tmop_pa_p3_c0.cpp
+7 -7
View File
@@ -41,7 +41,7 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_3D,
const auto C0 = const_c0 ?
Reshape(c0_.Read(), 1, 1, 1, 1) :
Reshape(c0_.Read(), Q1D, Q1D, Q1D, NE);
const auto LD = Reshape(lim_dist.Read(), D1D, D1D, D1D, DIM, NE);
const auto LD = Reshape(lim_dist.Read(), D1D, D1D, D1D, NE);
const auto J = Reshape(j_.Read(), DIM, DIM, Q1D, Q1D, Q1D, NE);
const auto b = Reshape(b_.Read(), Q1D, D1D);
const auto W = Reshape(w_.Read(), Q1D, Q1D, Q1D);
@@ -59,10 +59,10 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_3D,
MFEM_SHARED double B[MQ1*MD1];
MFEM_SHARED double DDD[3][MD1*MD1*MD1];
MFEM_SHARED double DDQ[3][MD1*MD1*MQ1];
MFEM_SHARED double DQQ[3][MD1*MQ1*MQ1];
MFEM_SHARED double QQQ[3][MQ1*MQ1*MQ1];
MFEM_SHARED double DDD[MD1*MD1*MD1];
MFEM_SHARED double DDQ[MD1*MD1*MQ1];
MFEM_SHARED double DQQ[MD1*MQ1*MQ1];
MFEM_SHARED double QQQ[MQ1*MQ1*MQ1];
MFEM_SHARED double DDD0[3][MD1*MD1*MD1];
MFEM_SHARED double DDQ0[3][MD1*MD1*MQ1];
@@ -102,7 +102,7 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_3D,
const double detJtr = kernels::Det<3>(Jtr);
const double weight = W(qx,qy,qz) * detJtr;
double D[3], p0[3], p1[3];
double D, p0[3], p1[3];
const double coeff0 = const_c0 ? C0(0,0,0,0) : C0(qx,qy,qz,e);
kernels::PullEval<MQ1>(qx,qy,qz,QQQ,D);
kernels::PullEval<MQ1>(qx,qy,qz,QQQ0,p0);
@@ -113,7 +113,7 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_3D,
// subtract(1.0 / (dist * dist), x, x0, d1);
// z = a * (x - y)
// grad = a * (x - x0)
const double dist = D[0]; // GetValues, default comp set to 0
const double dist = D; // GetValues, default comp set to 0
const double a = 1.0 / (dist * dist);
const double w = weight * lim_normal * coeff0;
kernels::Subtract<3>(w*a, p1, p0, d1);
fem/tmop_pa_w2_c0.cpp
+6 -6
View File
@@ -44,7 +44,7 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_2D,
const auto C0 = const_c0 ?
Reshape(c0_.Read(), 1, 1, 1) :
Reshape(c0_.Read(), Q1D, Q1D, NE);
const auto LD = Reshape(lim_dist.Read(), D1D, D1D, DIM, NE);
const auto LD = Reshape(lim_dist.Read(), D1D, D1D, NE);
const auto J = Reshape(j_.Read(), DIM, DIM, Q1D, Q1D, NE);
const auto b = Reshape(b_.Read(), Q1D, D1D);
const auto W = Reshape(w_.Read(), Q1D, Q1D);
@@ -63,9 +63,9 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_2D,
MFEM_SHARED double B[MQ1*MD1];
MFEM_SHARED double XY[2][NBZ][MD1*MD1];
MFEM_SHARED double DQ[2][NBZ][MD1*MQ1];
MFEM_SHARED double QQ[2][NBZ][MQ1*MQ1];
MFEM_SHARED double XY[NBZ][MD1*MD1];
MFEM_SHARED double DQ[NBZ][MD1*MQ1];
MFEM_SHARED double QQ[NBZ][MQ1*MQ1];
MFEM_SHARED double XY0[2][NBZ][MD1*MD1];
MFEM_SHARED double DQ0[2][NBZ][MD1*MQ1];
@@ -94,7 +94,7 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_2D,
{
MFEM_FOREACH_THREAD(qx,x,Q1D)
{
double ld[2], p0[2], p1[2];
double ld, p0[2], p1[2];
const double *Jtr = &J(0,0,qx,qy,e);
const double detJtr = kernels::Det<2>(Jtr);
const double weight = W(qx,qy) * detJtr;
@@ -102,7 +102,7 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_2D,
kernels::PullEval<MQ1,NBZ>(qx,qy,QQ,ld);
kernels::PullEval<MQ1,NBZ>(qx,qy,QQ0,p0);
kernels::PullEval<MQ1,NBZ>(qx,qy,QQ1,p1);
const double dist = ld[0]; // GetValues, default comp set to 0
const double dist = ld; // GetValues, default comp set to 0
const double id2 = 0.5 / (dist*dist);
const double dsq = kernels::DistanceSquared<2>(p1,p0) * id2;
E(qx,qy,e) = weight * lim_normal * dsq * coeff0;
fem/tmop_pa_w3_c0.cpp
+7 -7
View File
@@ -42,7 +42,7 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_3D,
const auto C0 = const_c0 ?
Reshape(c0_.Read(), 1, 1, 1, 1) :
Reshape(c0_.Read(), Q1D, Q1D, Q1D, NE);
const auto LD = Reshape(lim_dist.Read(), D1D, D1D, D1D, DIM, NE);
const auto LD = Reshape(lim_dist.Read(), D1D, D1D, D1D, NE);
const auto J = Reshape(j_.Read(), DIM, DIM, Q1D, Q1D, Q1D, NE);
const auto b = Reshape(b_.Read(), Q1D, D1D);
const auto W = Reshape(w_.Read(), Q1D, Q1D, Q1D);
@@ -60,10 +60,10 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_3D,
MFEM_SHARED double B[MQ1*MD1];
MFEM_SHARED double DDD[3][MD1*MD1*MD1];
MFEM_SHARED double DDQ[3][MD1*MD1*MQ1];
MFEM_SHARED double DQQ[3][MD1*MQ1*MQ1];
MFEM_SHARED double QQQ[3][MQ1*MQ1*MQ1];
MFEM_SHARED double DDD[MD1*MD1*MD1];
MFEM_SHARED double DDQ[MD1*MD1*MQ1];
MFEM_SHARED double DQQ[MD1*MQ1*MQ1];
MFEM_SHARED double QQQ[MQ1*MQ1*MQ1];
MFEM_SHARED double DDD0[3][MD1*MD1*MD1];
MFEM_SHARED double DDQ0[3][MD1*MD1*MQ1];
@@ -99,7 +99,7 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_3D,
{
MFEM_FOREACH_THREAD(qx,x,Q1D)
{
double D[3], p0[3], p1[3];
double D, p0[3], p1[3];
const double *Jtr = &J(0,0,qx,qy,qz,e);
const double detJtr = kernels::Det<3>(Jtr);
const double weight = W(qx,qy,qz) * detJtr;
@@ -109,7 +109,7 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_3D,
kernels::PullEval<MQ1>(qx,qy,qz,QQQ0,p0);
kernels::PullEval<MQ1>(qx,qy,qz,QQQ1,p1);
const double dist = D[0]; // GetValues, default comp set to 0
const double dist = D; // GetValues, default comp set to 0
const double id2 = 0.5 / (dist*dist);
const double dsq = kernels::DistanceSquared<3>(p1,p0) * id2;