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

Deprecate Vector implicit conversion to double*

Browse Source
This commit is contained in:
Will Pazner
2022-09-06 13:05:38 -07:00
parent 502d2f5dc7
commit 5a455dd141
21 changed files with 250 additions and 230 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fem/fe/fe_rt.cpp
+32 -30
View File
@@ -742,9 +742,9 @@ RT_TriangleElement::RT_TriangleElement(const int p)
for (int k = 0; k < dof; k++)
{
const IntegrationPoint &ip = Nodes.IntPoint(k);
poly1d.CalcBasis(p, ip.x, shape_x);
poly1d.CalcBasis(p, ip.y, shape_y);
poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l);
poly1d.CalcBasis(p, ip.x, shape_x.HostWrite());
poly1d.CalcBasis(p, ip.y, shape_y.HostWrite());
poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l.HostWrite());
const double *n_k = nk + 2*dof2nk[k];
o = 0;
@@ -776,9 +776,9 @@ void RT_TriangleElement::CalcVShape(const IntegrationPoint &ip,
DenseMatrix u(dof, dim);
#endif
poly1d.CalcBasis(p, ip.x, shape_x);
poly1d.CalcBasis(p, ip.y, shape_y);
poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l);
poly1d.CalcBasis(p, ip.x, shape_x.HostWrite());
poly1d.CalcBasis(p, ip.y, shape_y.HostWrite());
poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l.HostWrite());
int o = 0;
for (int j = 0; j <= p; j++)
@@ -810,9 +810,10 @@ void RT_TriangleElement::CalcDivShape(const IntegrationPoint &ip,
Vector divu(dof);
#endif
poly1d.CalcBasis(p, ip.x, shape_x, dshape_x);
poly1d.CalcBasis(p, ip.y, shape_y, dshape_y);
poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l, dshape_l);
poly1d.CalcBasis(p, ip.x, shape_x.HostWrite(), dshape_x.HostWrite());
poly1d.CalcBasis(p, ip.y, shape_y.HostWrite(), dshape_y.HostWrite());
poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l.HostWrite(),
dshape_l.HostWrite());
int o = 0;
for (int j = 0; j <= p; j++)
@@ -914,10 +915,10 @@ RT_TetrahedronElement::RT_TetrahedronElement(const int p)
for (int m = 0; m < dof; m++)
{
const IntegrationPoint &ip = Nodes.IntPoint(m);
poly1d.CalcBasis(p, ip.x, shape_x);
poly1d.CalcBasis(p, ip.y, shape_y);
poly1d.CalcBasis(p, ip.z, shape_z);
poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l);
poly1d.CalcBasis(p, ip.x, shape_x.HostWrite());
poly1d.CalcBasis(p, ip.y, shape_y.HostWrite());
poly1d.CalcBasis(p, ip.z, shape_z.HostWrite());
poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l.HostWrite());
const double *nm = nk + 3*dof2nk[m];
o = 0;
@@ -953,10 +954,10 @@ void RT_TetrahedronElement::CalcVShape(const IntegrationPoint &ip,
DenseMatrix u(dof, dim);
#endif
poly1d.CalcBasis(p, ip.x, shape_x);
poly1d.CalcBasis(p, ip.y, shape_y);
poly1d.CalcBasis(p, ip.z, shape_z);
poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l);
poly1d.CalcBasis(p, ip.x, shape_x.HostWrite());
poly1d.CalcBasis(p, ip.y, shape_y.HostWrite());
poly1d.CalcBasis(p, ip.z, shape_z.HostWrite());
poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l.HostWrite());
int o = 0;
for (int k = 0; k <= p; k++)
@@ -990,10 +991,11 @@ void RT_TetrahedronElement::CalcDivShape(const IntegrationPoint &ip,
Vector divu(dof);
#endif
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);
poly1d.CalcBasis(p, ip.x, shape_x.HostWrite(), dshape_x.HostWrite());
poly1d.CalcBasis(p, ip.y, shape_y.HostWrite(), dshape_y.HostWrite());
poly1d.CalcBasis(p, ip.z, shape_z.HostWrite(), dshape_z.HostWrite());
poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l.HostWrite(),
dshape_l.HostWrite());
int o = 0;
for (int k = 0; k <= p; k++)
@@ -1400,7 +1402,7 @@ void RT_R1D_SegmentElement::Project(const FiniteElement &fe,
Trans.SetIntPoint(&ip);
// Transform RT face normals from reference to physical space
// vk = adj(J)^T nk
Trans.AdjugateJacobian().MultTranspose(n1, vk);
Trans.AdjugateJacobian().MultTranspose(n1.HostRead(), vk);
vk[1] = n3[1] * Trans.Weight();
vk[2] = n3[2] * Trans.Weight();
if (fe.GetMapType() == INTEGRAL)
@@ -1446,7 +1448,7 @@ void RT_R1D_SegmentElement::Project(const FiniteElement &fe,
Trans.SetIntPoint(&ip);
// Transform RT face normals from reference to physical space
// vk = adj(J)^T nk
Trans.AdjugateJacobian().MultTranspose(n1, vk);
Trans.AdjugateJacobian().MultTranspose(n1.HostRead(), vk);
// Compute fe basis functions in physical space
fe.CalcVShape(Trans, vshape);
// Project fe basis functions onto transformed face normals
@@ -1477,7 +1479,7 @@ void RT_R1D_SegmentElement::ProjectCurl(const FiniteElement &fe,
for (int k = 0; k < dof; k++)
{
fe.CalcCurlShape(Nodes.IntPoint(k), curl_shape);
curl_shape.Mult(nk_ptr + dof2nk[k] * 3, curl_k);
curl_shape.Mult(nk_ptr + dof2nk[k] * 3, curl_k.HostWrite());
for (int j = 0; j < curl_k.Size(); j++)
{
curl(k,j) = (fabs(curl_k(j)) < 1e-12) ? 0.0 : curl_k(j);
@@ -1581,7 +1583,7 @@ void RT_R2D_SegmentElement::LocalInterpolation(const VectorFiniteElement &cfe,
ip.Set3(vk);
cfe.CalcVShape(ip, vshape);
// xk = |J| J^{-t} n_k
adjJ.MultTranspose(n2, vk);
adjJ.MultTranspose(n2.HostRead(), vk);
// I_k = vshape_k.adj(J)^t.n_k, k=1,...,dof
for (int j = 0; j < vshape.Height(); j++)
{
@@ -1654,7 +1656,7 @@ RT_R2D_FiniteElement::LocalInterpolation(const VectorFiniteElement &cfe,
ip.Set3(vk);
cfe.CalcVShape(ip, vshape);
// xk = |J| J^{-t} n_k
adjJ.MultTranspose(n2, vk);
adjJ.MultTranspose(n2.HostRead(), vk);
// I_k = vshape_k.adj(J)^t.n_k, k=1,...,dof
for (int j = 0; j < vshape.Height(); j++)
{
@@ -1695,7 +1697,7 @@ void RT_R2D_FiniteElement::GetLocalRestriction(ElementTransformation &Trans,
if (Geometries.CheckPoint(geom_type, ip)) // do we need an epsilon here?
{
CalcVShape(ip, vshape);
J.MultTranspose(n2, pt_data);
J.MultTranspose(n2.HostRead(), pt_data);
pt /= weight;
for (int k = 0; k < dof; k++)
{
@@ -1764,7 +1766,7 @@ void RT_R2D_FiniteElement::Project(const FiniteElement &fe,
Trans.SetIntPoint(&ip);
// Transform RT face normals from reference to physical space
// vk = adj(J)^T nk
Trans.AdjugateJacobian().MultTranspose(n2, vk);
Trans.AdjugateJacobian().MultTranspose(n2.HostRead(), vk);
vk[2] = n3[2] * Trans.Weight();
if (fe.GetMapType() == INTEGRAL)
{
@@ -1809,7 +1811,7 @@ void RT_R2D_FiniteElement::Project(const FiniteElement &fe,
Trans.SetIntPoint(&ip);
// Transform RT face normals from reference to physical space
// vk = adj(J)^T nk
Trans.AdjugateJacobian().MultTranspose(n2, vk);
Trans.AdjugateJacobian().MultTranspose(n2.HostRead(), vk);
// Compute fe basis functions in physical space
fe.CalcVShape(Trans, vshape);
// Project fe basis functions onto transformed face normals
@@ -1842,7 +1844,7 @@ void RT_R2D_FiniteElement::ProjectCurl(const FiniteElement &fe,
for (int k = 0; k < dof; k++)
{
fe.CalcCurlShape(Nodes.IntPoint(k), curl_shape);
curl_shape.Mult(nk_ptr + dof2nk[k] * 3, curl_k);
curl_shape.Mult(nk_ptr + dof2nk[k] * 3, curl_k.HostWrite());
for (int j = 0; j < curl_k.Size(); j++)
{
curl(k,j) = (fabs(curl_k(j)) < 1e-12) ? 0.0 : curl_k(j);