Skip to content

Commit

Permalink
Add support for p >= 2 Nedelec spaces
Browse files Browse the repository at this point in the history
  • Loading branch information
@sebastiangrimberg
sebastiangrimberg committed Apr 25, 2023
1 parent dfa39ef commit 19a4503
Show file tree
Hide file tree
Showing 2 changed files with 85 additions and 23 deletions.
103 changes: 83 additions & 20 deletions fem/ceed/interface/restriction.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -97,7 +97,7 @@ static void InitNativeRestr(const mfem::FiniteElementSpace &fes,

for (int i = 0; i < nelem; i++)
{
// TODO: Implement DofTransformation support
// DofTransformation support uses InitNativeRestrWithIndices
for (int j = 0; j < P; j++)
{
const int sgid = el_map[j + P * i]; // signed
Expand Down Expand Up @@ -149,7 +149,7 @@ static void InitLexicoRestrWithIndices(const mfem::FiniteElementSpace &fes,
{
// No need to handle DofTransformation for tensor-product elements
const int elem_index = indices[i];
DofTransformation *dof_trans;
mfem::DofTransformation *dof_trans;
if (use_bdr)
{
dof_trans = fes.GetBdrElementDofs(elem_index, dofs);
Expand Down Expand Up @@ -197,21 +197,31 @@ static void InitNativeRestrWithIndices(const mfem::FiniteElementSpace &fes,
Ceed ceed,
CeedElemRestriction *restr)
{
const mfem::FiniteElement *fe = use_bdr ? fes.GetBE(indices[0]) :
fes.GetFE(indices[0]);
const int i0 = indices ? indices[0] : 0;
const mfem::FiniteElement *fe = use_bdr ? fes.GetBE(i0) : fes.GetFE(i0);
const int P = fe->GetDof();
CeedInt compstride =
(fes.GetOrdering() == Ordering::byVDIM) ? 1 : fes.GetNDofs();
const int stride = (compstride == 1) ? fes.GetVDim() : 1;
mfem::Array<int> tp_el_dof(nelem * P), dofs;
mfem::Array<bool> tp_el_orients(nelem * P);
bool use_orients = false;
mfem::Array<bool> tp_el_orients;
mfem::Array<int> tp_el_curl_orients;
mfem::Vector el_trans_j;
mfem::DofTransformation *dof_trans = use_bdr ? fes.GetBdrElementDofs(i0, dofs) :
fes.GetElementDofs(i0, dofs);
if (!dof_trans)
{
tp_el_orients.SetSize(nelem * P);
}
else
{
tp_el_curl_orients.SetSize(nelem * P * 3, 0.0);
el_trans_j.SetSize(P);
}

for (int i = 0; i < nelem; i++)
{
// TODO: Implement DofTransformation support
const int elem_index = indices[i];
DofTransformation *dof_trans;
const int elem_index = indices ? indices[i] : i;
if (use_bdr)
{
dof_trans = fes.GetBdrElementDofs(elem_index, dofs);
Expand All @@ -220,19 +230,61 @@ static void InitNativeRestrWithIndices(const mfem::FiniteElementSpace &fes,
{
dof_trans = fes.GetElementDofs(elem_index, dofs);
}
MFEM_VERIFY(!dof_trans || fes.GetMaxElementOrder() == 1,
"DofTransformation support for CeedElemRestriction does not exist yet.");
for (int j = 0; j < P; j++)
if (!dof_trans)
{
const int sgid = dofs[j]; // signed
const int gid = (sgid >= 0) ? sgid : -1 - sgid;
tp_el_dof[j + P * i] = stride * gid;
tp_el_orients[j + P * i] = (sgid < 0);
use_orients = use_orients || tp_el_orients[j + P * i];
for (int j = 0; j < P; j++)
{
const int sgid = dofs[j]; // signed
const int gid = (sgid >= 0) ? sgid : -1 - sgid;
tp_el_dof[j + P * i] = stride * gid;
tp_el_orients[j + P * i] = (sgid < 0);
}
}
else
{
for (int j = 0; j < P; j++)
{
const int sgid = dofs[j]; // signed
const int gid = (sgid >= 0) ? sgid : -1 - sgid;
tp_el_dof[j + P * i] = stride * gid;

// Fill column j of element tridiagonal matrix tp_el_curl_orients
el_trans_j = 0.0;
el_trans_j(j) = 1.0;
dof_trans->InvTransformPrimal(el_trans_j);
el_trans_j *= (sgid < 0) ? -1.0 : 1.0;
tp_el_curl_orients[3 * (j + 0 + P * i) + 1] = el_trans_j(j + 0);
if (j > 0)
{
tp_el_curl_orients[3 * (j - 1 + P * i) + 2] = el_trans_j(j - 1);
}
if (j < P - 1)
{
tp_el_curl_orients[3 * (j + 1 + P * i) + 0] = el_trans_j(j + 1);
}
#ifdef MFEM_DEBUG
int nnz = 0;
for (int k = 0; k < P; k++)
{
if (k < j - 1 && k > j + 1 && el_trans_j(k) != 0.0) { nnz++; }
}
MFEM_ASSERT(nnz == 0,
"Element transformation matrix is not tridiagonal at column "
<< j << " (nnz = " << nnz << ")!");
#endif
}
}
}

if (use_orients)
if (tp_el_curl_orients.Size())
{
CeedElemRestrictionCreateCurlOriented(ceed, nelem, P, fes.GetVDim(),
compstride, fes.GetVDim() * fes.GetNDofs(),
CEED_MEM_HOST, CEED_COPY_VALUES,
tp_el_dof.GetData(), tp_el_curl_orients.GetData(),
restr);
}
else if (tp_el_orients.Size())
{
CeedElemRestrictionCreateOriented(ceed, nelem, P, fes.GetVDim(),
compstride, fes.GetVDim() * fes.GetNDofs(),
Expand All @@ -254,10 +306,13 @@ static void InitRestrictionImpl(const mfem::FiniteElementSpace &fes,
Ceed ceed,
CeedElemRestriction *restr)
{
const mfem::FiniteElement *fe = use_bdr ? fes.GetBE(0): fes.GetFE(0);
const mfem::FiniteElement *fe = use_bdr ? fes.GetBE(0) : fes.GetFE(0);
const mfem::TensorBasisElement *tfe =
dynamic_cast<const mfem::TensorBasisElement *>(fe);
const bool vector = fe->GetRangeType() == mfem::FiniteElement::VECTOR;
mfem::Array<int> dofs;
mfem::DofTransformation *dof_trans = use_bdr ? fes.GetBdrElementDofs(0, dofs) :
fes.GetElementDofs(0, dofs);
if (tfe && tfe->GetDofMap().Size() > 0 && !vector)
{
// Lexicographic ordering using dof_map
Expand All @@ -266,7 +321,15 @@ static void InitRestrictionImpl(const mfem::FiniteElementSpace &fes,
else
{
// Native ordering
InitNativeRestr(fes, use_bdr, ceed, restr);
if (!dof_trans)
{
InitNativeRestr(fes, use_bdr, ceed, restr);
}
else
{
InitNativeRestrWithIndices(fes, use_bdr, use_bdr ? fes.GetNBE() : fes.GetNE(),
nullptr, ceed, restr);
}
}
}

Expand Down
5 changes: 2 additions & 3 deletions tests/unit/ceed/test_ceed.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -985,14 +985,13 @@ TEST_CASE("CEED p-adaptivity", "[CEED]")
TEST_CASE("CEED vector and matrix coefficients and vector FE operators",
"[CEED], [VectorFE]")
{

auto assembly = GENERATE(AssemblyLevel::PARTIAL,AssemblyLevel::NONE);
auto coeff_type = GENERATE(CeedCoeffType::Const,CeedCoeffType::Quad,
CeedCoeffType::VecConst,CeedCoeffType::VecQuad,
CeedCoeffType::MatConst,CeedCoeffType::MatQuad);
auto pb = GENERATE(Problem::MassDiffusion,Problem::VectorFEMassDivDiv,
Problem::VectorFEMassCurlCurl);
auto order = GENERATE(1); // TODO p=2 curl-curl on tet mesh is not supported
auto order = GENERATE(1,3);
auto bdr_integ = GENERATE(false,true);
auto mesh = GENERATE("../../data/inline-quad.mesh",
"../../data/inline-hex.mesh",
Expand All @@ -1015,7 +1014,7 @@ TEST_CASE("CEED mixed integrators",
CeedCoeffType::VecConst,CeedCoeffType::VecQuad,
CeedCoeffType::MatConst,CeedCoeffType::MatQuad);
auto pb = GENERATE(Problem::MixedVectorGradient,Problem::MixedVectorCurl);
auto order = GENERATE(1); // TODO p=2 curl-curl on tet mesh is not supported
auto order = GENERATE(2);
auto bdr_integ = GENERATE(false,true);
auto mesh = GENERATE("../../data/inline-quad.mesh",
"../../data/inline-hex.mesh",
Expand Down

0 comments on commit 19a4503

Please sign in to comment.