From 509ad616fbe075c356fb56223624efce4244cd30 Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Thu, 11 Jan 2018 08:51:41 +0100
Subject: [PATCH] [boxLocalAssembler] Put treatment of Dirichlet BCs in
separate method
---
dumux/assembly/boxlocalassembler.hh | 74 ++++++++++++++++-------------
1 file changed, 41 insertions(+), 33 deletions(-)
diff --git a/dumux/assembly/boxlocalassembler.hh b/dumux/assembly/boxlocalassembler.hh
index e407740aa4..62e432f123 100644
--- a/dumux/assembly/boxlocalassembler.hh
+++ b/dumux/assembly/boxlocalassembler.hh
@@ -80,39 +80,7 @@ public:
for (const auto& scv : scvs(this->fvGeometry()))
res[scv.dofIndex()] += residual[scv.indexInElement()];
- // enforce Dirichlet boundaries by overwriting partial derivatives with 1 or 0
- // and set the residual to (privar - dirichletvalue)
- // TODO: put this in separate method!
- if (this->elemBcTypes().hasDirichlet())
- {
- for (const auto& scvI : scvs(this->fvGeometry()))
- {
- const auto bcTypes = this->elemBcTypes()[scvI.indexInElement()];
- if (bcTypes.hasDirichlet())
- {
- const auto dirichletValues = this->problem().dirichlet(this->element(), scvI);
-
- // set the dirichlet conditions in residual and jacobian
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
- {
- if (bcTypes.isDirichlet(eqIdx))
- {
- const auto pvIdx = bcTypes.eqToDirichletIndex(eqIdx);
- assert(0 <= pvIdx && pvIdx < numEq);
-
- const auto& priVars = this->curElemVolVars()[scvI].priVars();
- res[scvI.dofIndex()][eqIdx] = priVars[pvIdx] - dirichletValues[pvIdx];
- for (const auto& scvJ : scvs(this->fvGeometry()))
- {
- jac[scvI.dofIndex()][scvJ.dofIndex()][eqIdx] = 0.0;
- if (scvI.indexInElement() == scvJ.indexInElement())
- jac[scvI.dofIndex()][scvI.dofIndex()][eqIdx][pvIdx] = 1.0;
- }
- }
- }
- }
- }
- }
+ this->asImp_().evalDirichletBoundaries(jac, res);
}
/*!
@@ -123,6 +91,7 @@ public:
{
this->asImp_().bindLocalViews();
this->asImp_().assembleJacobianAndResidualImpl(jac, gridVariables); // forward to the internal implementation
+ // TODO: Dirichlet required here?
}
/*!
@@ -135,6 +104,8 @@ public:
for (const auto& scv : scvs(this->fvGeometry()))
res[scv.dofIndex()] += residual[scv.indexInElement()];
+
+ // TODO: Dirichlet required here?
}
@@ -147,6 +118,43 @@ public:
{
return this->evalLocalStorageResidual();
}
+
+ void evalDirichletBoundaries(JacobianMatrix& jac, SolutionVector& res)
+ {
+ // enforce Dirichlet boundaries by overwriting partial derivatives with 1 or 0
+ // and set the residual to (privar - dirichletvalue)
+ if (this->elemBcTypes().hasDirichlet())
+ {
+ for (const auto& scvI : scvs(this->fvGeometry()))
+ {
+ const auto bcTypes = this->elemBcTypes()[scvI.indexInElement()];
+ if (bcTypes.hasDirichlet())
+ {
+ const auto dirichletValues = this->problem().dirichlet(this->element(), scvI);
+
+ // set the dirichlet conditions in residual and jacobian
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ {
+ if (bcTypes.isDirichlet(eqIdx))
+ {
+ const auto pvIdx = bcTypes.eqToDirichletIndex(eqIdx);
+ assert(0 <= pvIdx && pvIdx < numEq);
+
+ const auto& priVars = this->curElemVolVars()[scvI].priVars();
+ res[scvI.dofIndex()][eqIdx] = priVars[pvIdx] - dirichletValues[pvIdx];
+ for (const auto& scvJ : scvs(this->fvGeometry()))
+ {
+ jac[scvI.dofIndex()][scvJ.dofIndex()][eqIdx] = 0.0;
+ if (scvI.indexInElement() == scvJ.indexInElement())
+ jac[scvI.dofIndex()][scvI.dofIndex()][eqIdx][pvIdx] = 1.0;
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
};
/*!
--
GitLab