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

implemented the LEXICOGRAPHIC_FULL idea

Browse Source
This commit is contained in:
Victor DeCaria
2023-12-04 16:55:50 -07:00
parent 4a38577e3f
commit 336236eae8
3 changed files with 162 additions and 23 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fem/fe/fe_base.cpp
+143 -12
View File
@@ -367,11 +367,10 @@ const DofToQuad &FiniteElement::GetDofToQuad(const IntegrationRule &ir,
{
MFEM_VERIFY(mode == DofToQuad::FULL, "invalid mode requested");
ElementDofOrdering ordering = ElementDofOrdering::NATIVE;
for (int i = 0; i < dof2quad_array.Size(); i++)
{
const DofToQuad &d2q = *dof2quad_array[i];
if (d2q.IntRule == &ir && d2q.mode == mode && d2q.ordering == ordering) { return d2q; }
if (d2q.IntRule == &ir && d2q.mode == mode) { return d2q; }
}
#ifdef MFEM_THREAD_SAFE
@@ -644,6 +643,118 @@ void ScalarFiniteElement::ScalarLocalL2Restriction(
}
}
void NodalFiniteElement::CreateLexicographicFullMap(const IntegrationRule &ir)
const
{
// Get the FULL version of the map.
auto &d2q = GetDofToQuad(ir, DofToQuad::FULL);
//Undo the native ordering which is the default in GetDofToQuad for FULL mode.
auto *d2q_new = new DofToQuad(d2q);
d2q_new->mode = DofToQuad::LEXICOGRAPHIC_FULL;
const int nqpt = ir.GetNPoints();
if (range_type == SCALAR)
{
for (int i = 0; i < nqpt; i++)
{
for (int j = 0; j < dof; j++)
{
d2q_new->B[i+nqpt*j] = d2q_new->Bt[j+dof*i] = d2q.B[i+nqpt*lex_ordering[j]];
}
}
}
else if (range_type == VECTOR)
{
for (int i = 0; i < nqpt; i++)
{
for (int d = 0; d < dim; d++)
{
for (int j = 0; j < dof; j++)
{
d2q_new->B[i+nqpt*(d+dim*j)] = d2q_new->Bt[j+dof*(i+nqpt*d)] = d2q.B[i+nqpt*
(d+dim*lex_ordering[j])];
}
}
}
}
else
{
// Skip B and Bt for unknown range type
}
switch (deriv_type)
{
case GRAD:
{
for (int i = 0; i < nqpt; i++)
{
for (int d = 0; d < dim; d++)
{
for (int j = 0; j < dof; j++)
{
d2q_new->G[i+nqpt*(d+dim*j)] = d2q_new->Gt[j+dof*(i+nqpt*d)] = d2q.G[i+nqpt*
(d+dim*lex_ordering[j])];
}
}
}
break;
}
case DIV:
{
for (int i = 0; i < nqpt; i++)
{
for (int j = 0; j < dof; j++)
{
d2q_new->G[i+nqpt*j] = d2q_new->Gt[j+dof*i] = d2q.G[i+nqpt*lex_ordering[j]];
}
}
break;
}
case CURL:
{
for (int i = 0; i < nqpt; i++)
{
for (int d = 0; d < cdim; d++)
{
for (int j = 0; j < dof; j++)
{
d2q_new->G[i+nqpt*(d+cdim*j)] = d2q_new->Gt[j+dof*(i+nqpt*d)] = d2q.G[i+nqpt*
(d+cdim*lex_ordering[j])];
}
}
}
break;
}
case NONE:
default:
// Skip G and Gt for unknown derivative type
break;
}
dof2quad_array.Append(d2q_new);
}
const DofToQuad &NodalFiniteElement::GetDofToQuad(const IntegrationRule &ir,
DofToQuad::Mode mode) const
{
MFEM_VERIFY(mode == DofToQuad::FULL ||
mode == DofToQuad::LEXICOGRAPHIC_FULL, "invalid mode requested");
//Should make this loop a function of FiniteElement
for (int i = 0; i < dof2quad_array.Size(); i++)
{
const DofToQuad &d2q = *dof2quad_array[i];
if (d2q.IntRule == &ir && d2q.mode == mode) { return d2q; }
}
if (mode == DofToQuad::FULL)
{
return FiniteElement::GetDofToQuad(ir, mode);
}
else
{
CreateLexicographicFullMap(ir);
return NodalFiniteElement::GetDofToQuad(ir, mode);
}
}
void NodalFiniteElement::ProjectCurl_2D(
const FiniteElement &fe, ElementTransformation &Trans,
DenseMatrix &curl) const
@@ -2532,29 +2643,49 @@ void NodalTensorFiniteElement::SetMapType(const int map_type)
const DofToQuad &NodalTensorFiniteElement::GetDofToQuad(
const IntegrationRule &ir,
DofToQuad::Mode mode,
const ElementDofOrdering ordering) const
DofToQuad::Mode mode) const
{
MFEM_VERIFY(!(mode == DofToQuad::Mode::TENSOR &&
ordering == ElementDofOrdering::NATIVE),
"Invalide combination of DofToQuad::Mode and ElementDofOrdering.");
MFEM_VERIFY(mode == DofToQuad::FULL ||
mode == DofToQuad::TENSOR ||
mode == DofToQuad::LEXICOGRAPHIC_FULL, "invalid mode requested");
//Should make this loop a function of FiniteElement
for (int i = 0; i < dof2quad_array.Size(); i++)
{
const DofToQuad &d2q = *dof2quad_array[i];
if (d2q.IntRule == &ir && d2q.mode == mode && d2q.ordering == ordering) { return d2q; }
if (d2q.IntRule == &ir && d2q.mode == mode) { return d2q; }
}
if (mode != DofToQuad::TENSOR)
{
return NodalFiniteElement::GetDofToQuad(ir, mode);
}
else
{
return GetTensorDofToQuad(*this, ir, mode, basis1d, true, dof2quad_array);
}
}
const DofToQuad &NodalTensorFiniteElement::GetDofToQuad(
const IntegrationRule &ir,
DofToQuad::Mode mode,
const ElementDofOrdering ordering) const
{
for (int i = 0; i < dof2quad_array.Size(); i++)
{
const DofToQuad &d2q = *dof2quad_array[i];
if (d2q.IntRule == &ir && d2q.mode == mode) { return d2q; }
}
//First create the DofToQuad map for ElementDofOrdering::NATIVE.
//Either return d2q, or create the FULL and LEXICOGRAPHIC map and return.
auto &d2q = GetDofToQuad(ir, mode);
if (mode == DofToQuad::Mode::FULL &&
ordering == ElementDofOrdering::LEXICOGRAPHIC)
if (mode == DofToQuad::LEXICOGRAPHIC_FULL)
{
//Undo the native ordering which is the default in GetDofToQuad for FULL mode.
auto *d2q_new = new DofToQuad(d2q);
d2q_new->ordering = ElementDofOrdering::LEXICOGRAPHIC;
d2q_new->mode = DofToQuad::LEXICOGRAPHIC_FULL;
const int nqpt = ir.GetNPoints();
if (range_type == SCALAR)
{
for (int i = 0; i < nqpt; i++)