Skip to content

Commit 2ee2569

Browse files
bsumirakbsumirak
committed
found an error in ConvectionDiffusion exports: with hanging nodes, the derivatives of the exports were not correctly assembled
1 parent 194a2ad commit 2ee2569

File tree

2 files changed

+79
-22
lines changed

2 files changed

+79
-22
lines changed

‎fv/convection_diffusion_fv.cpp‎

Lines changed: 32 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -908,8 +908,8 @@ ex_value(number vValue[],
908908
{
909909
// request for trial space
910910
try{
911-
const LocalShapeFunctionSet<dim>& rTrialSpace
912-
= LocalFiniteElementProvider::get<dim>(roid, m_lfeID);
911+
const LocalShapeFunctionSet<TFVGeom::dim>& rTrialSpace
912+
= LocalFiniteElementProvider::get<TFVGeom::dim>(roid, m_lfeID);
913913

914914
// storage for shape function at ip
915915
std::vector<number> vShape(rTrialSpace.num_sh());
@@ -993,8 +993,8 @@ ex_grad(MathVector<dim> vValue[],
993993
{
994994
// request for trial space
995995
try{
996-
const LocalShapeFunctionSet<dim>& rTrialSpace
997-
= LocalFiniteElementProvider::get<dim>(roid, m_lfeID);
996+
const LocalShapeFunctionSet<TFVGeom::dim>& rTrialSpace
997+
= LocalFiniteElementProvider::get<TFVGeom::dim>(roid, m_lfeID);
998998

999999
// storage for shape function at ip
10001000
const size_t numSH = rTrialSpace.num_sh();
@@ -1054,6 +1054,19 @@ register_all_funcs(const LFEID& lfeID, const int quadOrder)
10541054
const int order = lfeID.order();
10551055
if(quadOrder == order+1 && lfeID.type() == LFEID::LAGRANGE)
10561056
{
1057+
// RegularEdge
1058+
switch(order)
1059+
{
1060+
case 1: {typedef FVGeometry<1, RegularEdge, dim> FVGeom;
1061+
register_func<RegularEdge, FVGeom >(); break;}
1062+
case 2: {typedef FVGeometry<2, RegularEdge, dim> FVGeom;
1063+
register_func<RegularEdge, FVGeom >(); break;}
1064+
case 3: {typedef FVGeometry<3, RegularEdge, dim> FVGeom;
1065+
register_func<RegularEdge, FVGeom >(); break;}
1066+
default: {typedef DimFVGeometry<dim, 1> FVGeom;
1067+
register_func<RegularEdge, FVGeom >(); break;}
1068+
}
1069+
10571070
// Triangle
10581071
switch(order)
10591072
{
@@ -1081,6 +1094,8 @@ register_all_funcs(const LFEID& lfeID, const int quadOrder)
10811094
}
10821095
else
10831096
{
1097+
register_func<RegularEdge, DimFVGeometry<dim, 1> >();
1098+
10841099
typedef DimFVGeometry<dim> FVGeom;
10851100
register_func<Triangle, FVGeom >();
10861101
register_func<Quadrilateral, FVGeom >();
@@ -1096,6 +1111,19 @@ register_all_funcs(const LFEID& lfeID, const int quadOrder)
10961111
const int order = lfeID.order();
10971112
if(quadOrder == order+1 && lfeID.type() == LFEID::LAGRANGE)
10981113
{
1114+
// RegularEdge
1115+
switch(order)
1116+
{
1117+
case 1: {typedef FVGeometry<1, RegularEdge, dim> FVGeom;
1118+
register_func<RegularEdge, FVGeom >(); break;}
1119+
case 2: {typedef FVGeometry<2, RegularEdge, dim> FVGeom;
1120+
register_func<RegularEdge, FVGeom >(); break;}
1121+
case 3: {typedef FVGeometry<3, RegularEdge, dim> FVGeom;
1122+
register_func<RegularEdge, FVGeom >(); break;}
1123+
default: {typedef DimFVGeometry<dim, 1> FVGeom;
1124+
register_func<RegularEdge, FVGeom >(); break;}
1125+
}
1126+
10991127
// Tetrahedron
11001128
switch(order)
11011129
{

‎fv1/convection_diffusion_fv1.cpp‎

Lines changed: 47 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -1395,6 +1395,7 @@ ex_value(number vValue[],
13951395
// number of shape functions
13961396
static const size_t numSH = ref_elem_type::numCorners;
13971397

1398+
13981399
// FV1 SCVF ip
13991400
if(vLocIP == geo.scvf_local_ips())
14001401
{
@@ -1411,22 +1412,41 @@ ex_value(number vValue[],
14111412

14121413
// compute derivative w.r.t. to unknowns iff needed
14131414
if(bDeriv)
1415+
{
14141416
for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
14151417
vvvDeriv[ip][_C_][sh] = scvf.shape(sh);
1418+
1419+
// do not forget that number of DoFs (== vvvDeriv[ip][_C_])
1420+
// might be > scvf.num_sh() in case of hanging nodes!
1421+
size_t ndof = vvvDeriv[ip][_C_].size();
1422+
for (size_t sh = scvf.num_sh(); sh < ndof; ++sh)
1423+
vvvDeriv[ip][_C_][sh] = 0.0;
1424+
}
14161425
}
14171426
}
14181427
// FV1 SCV ip
14191428
else if(vLocIP == geo.scv_local_ips())
14201429
{
1421-
// solution at ip
1422-
for(size_t sh = 0; sh < numSH; ++sh)
1423-
vValue[sh] = u(_C_, sh);
1430+
// Loop Sub Control Volumes (SCV)
1431+
for(size_t ip = 0; ip < geo.num_scv(); ++ip)
1432+
{
1433+
// Get current SCV
1434+
const typename TFVGeom::SCV& scv = geo.scv(ip);
14241435

1425-
// set derivatives if needed
1426-
if(bDeriv)
1427-
for(size_t sh = 0; sh < numSH; ++sh)
1428-
for(size_t sh2 = 0; sh2 < numSH; ++sh2)
1429-
vvvDeriv[sh][_C_][sh2] = (sh==sh2) ? 1.0 : 0.0;
1436+
// get corner of SCV
1437+
const size_t co = scv.node_id();
1438+
1439+
// solution at ip
1440+
vValue[ip] = u(_C_, co);
1441+
1442+
// set derivatives if needed
1443+
if(bDeriv)
1444+
{
1445+
size_t ndof = vvvDeriv[ip][_C_].size();
1446+
for(size_t sh = 0; sh < ndof; ++sh)
1447+
vvvDeriv[ip][_C_][sh] = (sh==co) ? 1.0 : 0.0;
1448+
}
1449+
}
14301450
}
14311451
// general case
14321452
else
@@ -1451,8 +1471,15 @@ ex_value(number vValue[],
14511471
// compute derivative w.r.t. to unknowns iff needed
14521472
// \todo: maybe store shapes directly in vvvDeriv
14531473
if(bDeriv)
1474+
{
14541475
for(size_t sh = 0; sh < numSH; ++sh)
14551476
vvvDeriv[ip][_C_][sh] = vShape[sh];
1477+
1478+
// beware of hanging nodes!
1479+
size_t ndof = vvvDeriv[ip][_C_].size();
1480+
for (size_t sh = numSH; sh < ndof; ++sh)
1481+
vvvDeriv[ip][_C_][sh] = 0.0;
1482+
}
14561483
}
14571484
}
14581485
}
@@ -1485,16 +1512,6 @@ ex_grad(MathVector<dim> vValue[],
14851512
// number of shape functions
14861513
static const size_t numSH = ref_elem_type::numCorners;
14871514

1488-
// reset the values for the derivatives
1489-
// this is necessary as vvvDeriv comes uninitialized
1490-
// and in the case of hanging nodes, vvvDeriv[ip][c].size() may be
1491-
// larger than geo.scvf(ip).num_sh(), thus some uninit'ed values
1492-
// are never changed (resulting in solver breakdowns, NaNs etc.)
1493-
if (bDeriv)
1494-
for (size_t ip = 0; ip < nip; ++ip)
1495-
for (size_t c = 0; c < vvvDeriv[ip].size(); ++c)
1496-
for (size_t sh = 0; sh < vvvDeriv[ip][c].size(); ++sh)
1497-
vvvDeriv[ip][c][sh] = 0.0;
14981515

14991516
// FV1 SCVF ip
15001517
if(vLocIP == geo.scvf_local_ips())
@@ -1514,6 +1531,11 @@ ex_grad(MathVector<dim> vValue[],
15141531
{
15151532
for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
15161533
vvvDeriv[ip][_C_][sh] = scvf.global_grad(sh);
1534+
1535+
// beware of hanging nodes!
1536+
size_t ndof = vvvDeriv[ip][_C_].size();
1537+
for (size_t sh = scvf.num_sh(); sh < ndof; ++sh)
1538+
vvvDeriv[ip][_C_][sh] = 0.0;
15171539
}
15181540
}
15191541
}
@@ -1548,8 +1570,15 @@ ex_grad(MathVector<dim> vValue[],
15481570

15491571
// compute derivative w.r.t. to unknowns iff needed
15501572
if(bDeriv)
1573+
{
15511574
for(size_t sh = 0; sh < numSH; ++sh)
15521575
MatVecMult(vvvDeriv[ip][_C_][sh], JTInv, vLocGrad[sh]);
1576+
1577+
// beware of hanging nodes!
1578+
size_t ndof = vvvDeriv[ip][_C_].size();
1579+
for (size_t sh = numSH; sh < ndof; ++sh)
1580+
vvvDeriv[ip][_C_][sh] = 0.0;
1581+
}
15531582
}
15541583
}
15551584
};

0 commit comments

Comments
 (0)