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

Merge remote-tracking branch 'origin/master' into pyramid-dev

Browse Source 
# Conflicts:
#	fem/fe/fe_base.cpp
#	fem/fe/fe_base.hpp
#	fem/gridfunc.cpp
This commit is contained in:
Stowell, Mark L
2024-04-08 12:09:46 -07:00
parent 1a2b8291f5 91710623f8
commit 4804824873
1120 changed files with 53619 additions and 24294 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fem/fe/fe_l2.cpp
+36 -36
View File
@@ -1,4 +1,4 @@
// Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
// LICENSE and NOTICE for details. LLNL-CODE-806117.
//
@@ -23,7 +23,7 @@ using namespace std;
L2_SegmentElement::L2_SegmentElement(const int p, const int btype)
: NodalTensorFiniteElement(1, p, VerifyOpen(btype), L2_DOF_MAP)
{
const double *op = poly1d.OpenPoints(p, btype);
const real_t *op = poly1d.OpenPoints(p, btype);
#ifndef MFEM_THREAD_SAFE
shape_x.SetSize(p + 1);
@@ -58,7 +58,7 @@ void L2_SegmentElement::CalcDShape(const IntegrationPoint &ip,
void L2_SegmentElement::ProjectDelta(int vertex, Vector &dofs) const
{
const int p = order;
const double *op = poly1d.OpenPoints(p, b_type);
const real_t *op = poly1d.OpenPoints(p, b_type);
switch (vertex)
{
@@ -82,7 +82,7 @@ void L2_SegmentElement::ProjectDelta(int vertex, Vector &dofs) const
L2_QuadrilateralElement::L2_QuadrilateralElement(const int p, const int btype)
: NodalTensorFiniteElement(2, p, VerifyOpen(btype), L2_DOF_MAP)
{
const double *op = poly1d.OpenPoints(p, b_type);
const real_t *op = poly1d.OpenPoints(p, b_type);
#ifndef MFEM_THREAD_SAFE
shape_x.SetSize(p + 1);
@@ -142,7 +142,7 @@ void L2_QuadrilateralElement::CalcDShape(const IntegrationPoint &ip,
void L2_QuadrilateralElement::ProjectDelta(int vertex, Vector &dofs) const
{
const int p = order;
const double *op = poly1d.OpenPoints(p, b_type);
const real_t *op = poly1d.OpenPoints(p, b_type);
#ifdef MFEM_THREAD_SAFE
Vector shape_x(p+1), shape_y(p+1);
@@ -200,16 +200,16 @@ void L2_QuadrilateralElement::ProjectDiv(const FiniteElement &fe,
div = 0.0;
const IntegrationRule &ir = IntRules.Get(geom_type, fe.GetOrder());
const double *gll_pts = poly1d.GetPoints(order+1, BasisType::GaussLobatto);
const real_t *gll_pts = poly1d.GetPoints(order+1, BasisType::GaussLobatto);
// Loop over subcells
for (int iy = 0; iy < order+1; ++iy)
{
const double hy = gll_pts[iy+1] - gll_pts[iy];
const real_t hy = gll_pts[iy+1] - gll_pts[iy];
for (int ix = 0; ix < order+1; ++ix)
{
const int i = ix + iy*(order+1);
const double hx = gll_pts[ix+1] - gll_pts[ix];
const real_t hx = gll_pts[ix+1] - gll_pts[ix];
// Loop over subcell quadrature points
for (int iq = 0; iq < ir.Size(); ++iq)
{
@@ -218,10 +218,10 @@ void L2_QuadrilateralElement::ProjectDiv(const FiniteElement &fe,
ip.y = gll_pts[iy] + hy*ip.y;
Trans.SetIntPoint(&ip);
fe.CalcDivShape(ip, div_shape);
double w = ip.weight;
real_t w = ip.weight;
if (map_type == VALUE)
{
const double detJ = Trans.Weight();
const real_t detJ = Trans.Weight();
w /= detJ;
}
else if (map_type == INTEGRAL)
@@ -230,7 +230,7 @@ void L2_QuadrilateralElement::ProjectDiv(const FiniteElement &fe,
}
for (int j = 0; j < fe_ndof; j++)
{
const double div_j = div_shape(j);
const real_t div_j = div_shape(j);
div(i,j) += w*div_j;
}
}
@@ -259,17 +259,17 @@ void L2_QuadrilateralElement::Project(Coefficient &coeff,
if (basis1d.IsIntegratedType())
{
const IntegrationRule &ir = IntRules.Get(geom_type, order);
const double *gll_pts = poly1d.GetPoints(order+1, BasisType::GaussLobatto);
const real_t *gll_pts = poly1d.GetPoints(order+1, BasisType::GaussLobatto);
dofs = 0.0;
// Loop over subcells
for (int iy = 0; iy < order+1; ++iy)
{
const double hy = gll_pts[iy+1] - gll_pts[iy];
const real_t hy = gll_pts[iy+1] - gll_pts[iy];
for (int ix = 0; ix < order+1; ++ix)
{
const int i = ix + iy*(order+1);
const double hx = gll_pts[ix+1] - gll_pts[ix];
const real_t hx = gll_pts[ix+1] - gll_pts[ix];
// Loop over subcell quadrature points
for (int iq = 0; iq < ir.Size(); ++iq)
{
@@ -277,8 +277,8 @@ void L2_QuadrilateralElement::Project(Coefficient &coeff,
ip.x = gll_pts[ix] + hx*ip.x;
ip.y = gll_pts[iy] + hy*ip.y;
Trans.SetIntPoint(&ip);
const double val = coeff.Eval(Trans, ip);
double w = ip.weight;
const real_t val = coeff.Eval(Trans, ip);
real_t w = ip.weight;
if (map_type == INTEGRAL)
{
w *= hx*hy*Trans.Weight();
@@ -298,7 +298,7 @@ void L2_QuadrilateralElement::Project(Coefficient &coeff,
L2_HexahedronElement::L2_HexahedronElement(const int p, const int btype)
: NodalTensorFiniteElement(3, p, VerifyOpen(btype), L2_DOF_MAP)
{
const double *op = poly1d.OpenPoints(p, btype);
const real_t *op = poly1d.OpenPoints(p, btype);
#ifndef MFEM_THREAD_SAFE
shape_x.SetSize(p + 1);
@@ -367,7 +367,7 @@ void L2_HexahedronElement::CalcDShape(const IntegrationPoint &ip,
void L2_HexahedronElement::ProjectDelta(int vertex, Vector &dofs) const
{
const int p = order;
const double *op = poly1d.OpenPoints(p, b_type);
const real_t *op = poly1d.OpenPoints(p, b_type);
#ifdef MFEM_THREAD_SAFE
Vector shape_x(p+1), shape_y(p+1);
@@ -461,19 +461,19 @@ void L2_HexahedronElement::ProjectDiv(const FiniteElement &fe,
div = 0.0;
const IntegrationRule &ir = IntRules.Get(geom_type, fe.GetOrder());
const double *gll_pts = poly1d.GetPoints(order+1, BasisType::GaussLobatto);
const real_t *gll_pts = poly1d.GetPoints(order+1, BasisType::GaussLobatto);
// Loop over subcells
for (int iz = 0; iz < order+1; ++iz)
{
const double hz = gll_pts[iz+1] - gll_pts[iz];
const real_t hz = gll_pts[iz+1] - gll_pts[iz];
for (int iy = 0; iy < order+1; ++iy)
{
const double hy = gll_pts[iy+1] - gll_pts[iy];
const real_t hy = gll_pts[iy+1] - gll_pts[iy];
for (int ix = 0; ix < order+1; ++ix)
{
const int i = ix + iy*(order+1) + iz*(order+1)*(order+1);
const double hx = gll_pts[ix+1] - gll_pts[ix];
const real_t hx = gll_pts[ix+1] - gll_pts[ix];
// Loop over subcell quadrature points
for (int iq = 0; iq < ir.Size(); ++iq)
{
@@ -483,10 +483,10 @@ void L2_HexahedronElement::ProjectDiv(const FiniteElement &fe,
ip.z = gll_pts[iz] + hz*ip.z;
Trans.SetIntPoint(&ip);
fe.CalcDivShape(ip, div_shape);
double w = ip.weight;
real_t w = ip.weight;
if (map_type == VALUE)
{
const double detJ = Trans.Weight();
const real_t detJ = Trans.Weight();
w /= detJ;
}
else if (map_type == INTEGRAL)
@@ -495,7 +495,7 @@ void L2_HexahedronElement::ProjectDiv(const FiniteElement &fe,
}
for (int j = 0; j < fe_ndof; j++)
{
const double div_j = div_shape(j);
const real_t div_j = div_shape(j);
div(i,j) += w*div_j;
}
}
@@ -525,19 +525,19 @@ void L2_HexahedronElement::Project(Coefficient &coeff,
if (basis1d.IsIntegratedType())
{
const IntegrationRule &ir = IntRules.Get(geom_type, order);
const double *gll_pts = poly1d.GetPoints(order+1, BasisType::GaussLobatto);
const real_t *gll_pts = poly1d.GetPoints(order+1, BasisType::GaussLobatto);
dofs = 0.0;
// Loop over subcells
for (int iz = 0; iz < order+1; ++iz)
{
const double hz = gll_pts[iz+1] - gll_pts[iz];
const real_t hz = gll_pts[iz+1] - gll_pts[iz];
for (int iy = 0; iy < order+1; ++iy)
{
const double hy = gll_pts[iy+1] - gll_pts[iy];
const real_t hy = gll_pts[iy+1] - gll_pts[iy];
for (int ix = 0; ix < order+1; ++ix)
{
const double hx = gll_pts[ix+1] - gll_pts[ix];
const real_t hx = gll_pts[ix+1] - gll_pts[ix];
const int i = ix + iy*(order+1) + iz*(order+1)*(order+1);
// Loop over subcell quadrature points
for (int iq = 0; iq < ir.Size(); ++iq)
@@ -547,11 +547,11 @@ void L2_HexahedronElement::Project(Coefficient &coeff,
ip.y = gll_pts[iy] + hy*ip.y;
ip.z = gll_pts[iz] + hz*ip.z;
Trans.SetIntPoint(&ip);
const double val = coeff.Eval(Trans, ip);
double w = ip.weight;
const real_t val = coeff.Eval(Trans, ip);
real_t w = ip.weight;
if (map_type == INTEGRAL)
{
const double detJ = Trans.Weight();
const real_t detJ = Trans.Weight();
w *= detJ*hx*hy*hz;
}
dofs[i] += val*w;
@@ -571,7 +571,7 @@ L2_TriangleElement::L2_TriangleElement(const int p, const int btype)
: NodalFiniteElement(2, Geometry::TRIANGLE, ((p + 1)*(p + 2))/2, p,
FunctionSpace::Pk)
{
const double *op = poly1d.OpenPoints(p, VerifyOpen(btype));
const real_t *op = poly1d.OpenPoints(p, VerifyOpen(btype));
#ifndef MFEM_THREAD_SAFE
shape_x.SetSize(p + 1);
@@ -589,7 +589,7 @@ L2_TriangleElement::L2_TriangleElement(const int p, const int btype)
for (int o = 0, j = 0; j <= p; j++)
for (int i = 0; i + j <= p; i++)
{
double w = op[i] + op[j] + op[p-i-j];
real_t w = op[i] + op[j] + op[p-i-j];
Nodes.IntPoint(o++).Set2(op[i]/w, op[j]/w);
}
@@ -696,7 +696,7 @@ L2_TetrahedronElement::L2_TetrahedronElement(const int p, const int btype)
: NodalFiniteElement(3, Geometry::TETRAHEDRON, ((p + 1)*(p + 2)*(p + 3))/6,
p, FunctionSpace::Pk)
{
const double *op = poly1d.OpenPoints(p, VerifyOpen(btype));
const real_t *op = poly1d.OpenPoints(p, VerifyOpen(btype));
#ifndef MFEM_THREAD_SAFE
shape_x.SetSize(p + 1);
@@ -717,7 +717,7 @@ L2_TetrahedronElement::L2_TetrahedronElement(const int p, const int btype)
for (int j = 0; j + k <= p; j++)
for (int i = 0; i + j + k <= p; i++)
{
double w = op[i] + op[j] + op[k] + op[p-i-j-k];
real_t w = op[i] + op[j] + op[k] + op[p-i-j-k];
Nodes.IntPoint(o++).Set3(op[i]/w, op[j]/w, op[k]/w);
}