Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

code for obtaining subcell fluxes with desired sparsity pattern #35

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

code for obtaining subcell fluxes with desired sparsity pattern #35

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
24 changes: 24 additions & 0 deletions remhos.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -854,6 +854,30 @@ int main(int argc, char *argv[])
s_max_glob = -numeric_limits<double>::infinity();
const int NE = pmesh.GetNE();

const FiniteElement *dummy = pfes.GetFE(0);
const int nd = dummy->GetDof();

DenseMatrix DistMat, MassMatLOR;
GetDistOps(dummy, DistMat, MassMatLOR);

Vector elem_vec(nd), elem_vec_dist(nd);
elem_vec.Randomize();
double aux = elem_vec.Sum() / nd;
elem_vec -= aux; // sum of vector entries must be zero

ApplyDistOp(DistMat, elem_vec, elem_vec_dist);

// Testing
for (int i = 0; i < nd; i++)
{
for (int j = 0; j < nd; j++)
{
double f_ij = GetAntiDiffFlux(MassMatLOR, elem_vec_dist, i, j);
cout << f_ij << endl;
}
cout << endl;
}

// Time-integration (loop over the time iterations, ti, with a time-step dt).
bool done = false;
for (int ti = 0; !done;)
Expand Down
155 changes: 155 additions & 0 deletions remhos_tools.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,161 @@ using namespace std;
namespace mfem
{

void ComputeLORMassMatRef(Geometry::Type gtype, bool UseDiagonalNbrs, DenseMatrix &RefMat)
{
switch (gtype)
{
case Geometry::SEGMENT:
{
RefMat = 1.0;
RefMat(0,0) = RefMat(1,1) = 2.0;
RefMat *= 1./6.;
break;
}
case Geometry::SQUARE:
{
RefMat = 2.0;
if (UseDiagonalNbrs)
{
RefMat(0,0) = RefMat(1,1) = RefMat(2,2) = RefMat(3,3) = 4.0;
RefMat(0,3) = RefMat(1,2) = RefMat(2,1) = RefMat(3,0) = 1.0;
}
else
{
RefMat(0,0) = RefMat(1,1) = RefMat(2,2) = RefMat(3,3) = 5.0;
RefMat(0,3) = RefMat(1,2) = RefMat(2,1) = RefMat(3,0) = 0.0;
}
RefMat *= 1.0/36.0;
break;
}
case Geometry::CUBE:
{
if (UseDiagonalNbrs)
{
RefMat = 2.0;
RefMat(0,0) = RefMat(1,1) = RefMat(2,2) = RefMat(3,3) = 8.0;
RefMat(4,4) = RefMat(5,5) = RefMat(6,6) = RefMat(7,7) = 8.0;
RefMat(0,7) = RefMat(1,6) = RefMat(2,5) = RefMat(3,4) = 1.0;
RefMat(4,3) = RefMat(5,2) = RefMat(6,1) = RefMat(7,0) = 1.0;
}
else
{
RefMat = 0.;
RefMat(0,0) = RefMat(1,1) = RefMat(2,2) = RefMat(3,3) = 15.0;
RefMat(4,4) = RefMat(5,5) = RefMat(6,6) = RefMat(7,7) = 15.0;
}

RefMat(0,1) = RefMat(0,2) = RefMat(1,0) = RefMat(1,3) = 4.0;
RefMat(2,0) = RefMat(2,3) = RefMat(3,1) = RefMat(3,2) = 4.0;
RefMat(4,5) = RefMat(4,6) = RefMat(5,4) = RefMat(5,7) = 4.0;
RefMat(6,4) = RefMat(6,7) = RefMat(7,5) = RefMat(7,6) = 4.0;
RefMat(0,4) = RefMat(1,5) = RefMat(2,6) = RefMat(3,7) = 4.0;
RefMat(4,0) = RefMat(5,1) = RefMat(6,2) = RefMat(7,3) = 4.0;

RefMat *= 1.0/216.0;
break;
}
default:
MFEM_ABORT("Unsupported Geometry.");
}
}

void GetDistOps(const FiniteElement *dummy, DenseMatrix &DistMat,
DenseMatrix &MassMatLOR)
{
const int p = dummy->GetOrder();
const int dim = dummy->GetDim();
const int nd = dummy->GetDof();
const Geometry::Type gtype = dummy->GetGeomType();
const int ndSub = int(pow(2, dim)); // number of closest neighbors
const int nSub = int(pow(p, dim)); // number of subcells

DenseMatrix Sub2Ind(nSub, ndSub);

for (int m = 0; m < nSub; m++)
{
for (int j = 0; j < ndSub; j++)
{
if (dim == 1) { Sub2Ind(m,j) = m + j; }
else if (dim == 2 && gtype == Geometry::SQUARE)
{
int aux = m + m/p;
switch (j)
{
case 0: Sub2Ind(m,j) = aux; break;
case 1: Sub2Ind(m,j) = aux + 1; break;
case 2: Sub2Ind(m,j) = aux + p+1; break;
case 3: Sub2Ind(m,j) = aux + p+2; break;
}
}
else if (dim == 3 && gtype == Geometry::CUBE)
{
int aux = m + m/p + (p+1)*(m/(p*p));
switch (j)
{
case 0: Sub2Ind(m,j) = aux; break;
case 1: Sub2Ind(m,j) = aux + 1; break;
case 2: Sub2Ind(m,j) = aux + p+1; break;
case 3: Sub2Ind(m,j) = aux + p+2; break;
case 4: Sub2Ind(m,j) = aux + (p+1)*(p+1); break;
case 5: Sub2Ind(m,j) = aux + (p+1)*(p+1)+1; break;
case 6: Sub2Ind(m,j) = aux + (p+1)*(p+1)+p+1; break;
case 7: Sub2Ind(m,j) = aux + (p+1)*(p+1)+p+2; break;
}
}
else
{
MFEM_ABORT("Simplices not required(?).");
}
}
}

MassMatLOR.SetSize(nd, nd);
MassMatLOR = 0.0;
DenseMatrix RefMat(ndSub, ndSub);
ComputeLORMassMatRef(gtype, false, RefMat);

for (int m = 0; m < nSub; m++)
{
for (int i = 0; i < ndSub; i++)
{
int I = Sub2Ind(m,i);
for (int j = 0; j < ndSub; j++)
{
int J = Sub2Ind(m,j);
MassMatLOR(I,J) += RefMat(i,j);
}
}
}

Vector MassMatLORLumped(nd);
DenseMatrix LowOrderRefinedMat(MassMatLOR);
LowOrderRefinedMat *= -1.0;
LowOrderRefinedMat.GetRowSums(MassMatLORLumped);

for (int i = 0; i < nd; i++)
{
LowOrderRefinedMat(i,i) -= MassMatLORLumped(i);
LowOrderRefinedMat(0,i) = 1.0;
}

DenseMatrixInverse InvLowOrderRefinedMat(&LowOrderRefinedMat);
InvLowOrderRefinedMat.Factor();
InvLowOrderRefinedMat.GetInverseMatrix(DistMat);
}

void ApplyDistOp(const DenseMatrix &DistMat, const Vector &elem_vec,
Vector &elem_vec_dist)
{
DistMat.Mult(elem_vec, elem_vec_dist);
}

double GetAntiDiffFlux(const DenseMatrix &MassMatLOR,
Vector &elem_vec_dist, int i, int j)
{
return MassMatLOR(i,j) * (elem_vec_dist(i) - elem_vec_dist(j));
}

SmoothnessIndicator::SmoothnessIndicator(int type_id,
ParMesh &subcell_mesh,
ParFiniteElementSpace &pfes_DG_,
Expand Down
9 changes: 9 additions & 0 deletions remhos_tools.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,15 @@
namespace mfem
{

void ComputeLORMassMatRef(Geometry::Type gtype, bool UseDiagonalNbrs,
DenseMatrix &RefMat);
void GetDistOps(const FiniteElement *dummy, DenseMatrix &DistMat,
DenseMatrix &MassMatLOR);
void ApplyDistOp(const DenseMatrix &DistMat, const Vector &elem_vec,
Vector &elem_vec_dist);
double GetAntiDiffFlux(const DenseMatrix &MassMatLOR,
Vector &elem_vec_dist, int i, int j);

int GetLocalFaceDofIndex(int dim, int loc_face_id, int face_orient,
int face_dof_id, int face_dof1D_cnt);

Expand Down