Commit c40b0e9a authored by Christian Engwer's avatar Christian Engwer
Browse files

[localized subtraction] fix several smaller bugs, as discussed with Alexander

parent c5bf64d4
......@@ -86,12 +86,12 @@ namespace duneuro
// sigma_
auto sigma = param.A(eg, Dune::ReferenceElements<DF, dim>::general(gt).position(0, 0));
//sigma_corr
auto sigma_corr = param.A(eg, Dune::ReferenceElements<DF, dim>::general(gt).position(0, 0));
sigma_corr -= param.get_sigma_infty();
//sigma_infty
auto sigma_infty = sigma;
auto sigma_infty = param.get_sigma_infty();
//sigma_corr
auto sigma_corr = sigma;
sigma_corr -= sigma_infty;
std::vector<RangeType> phi(lfsv.size());
std::vector<Dune::FieldMatrix<RF, 1, dim>> gradphi(lfsv.size());
......@@ -100,9 +100,9 @@ namespace duneuro
const auto& rule = Dune::QuadratureRules<DF, dim>::rule(gt, intorder);
for (const auto& qp : rule) {
BasisSwitch::gradient(FESwitch::basis(lfsv.finiteElement()), geometry, qp.position(),
gradphi);
gradphi);
//The volume term on the right hand side consits of three parts:
// The volume term on the right hand side consits of three parts:
// 1. part of the right hand side: (sigma u_infty) gradient_chi
auto sigma_u_infty = sigma;
......@@ -110,19 +110,19 @@ namespace duneuro
typename PROBLEMDATA::Traits::RangeType sigma_u_infty_grad_chi;
sigma_u_infty.mv(param.get_grad_chi(geometry.global(qp.position())),
sigma_u_infty_grad_chi);
// 2. part: sigma_corr grad_u_infty chi
typename PROBLEMDATA::Traits::RangeType sigma_corr_grad_u_infty_chi;
sigma_corr.mv(param.get_grad_u_infty(geometry.global(qp.position())), // see FieldMatrix, FieldVector
sigma_corr_grad_u_infty_chi);
sigma_corr_grad_u_infty_chi *= param.get_chi(qp.position());
// 3. part: sigma_infty grad_u_infty (1 - chi)
typename PROBLEMDATA::Traits::RangeType sigma_infty_grad_u_infty;
sigma_infty.mv(param.get_grad_u_infty(geometry.global(qp.position())), // see FieldMatrix, FieldVector
sigma_infty_grad_u_infty);
typename PROBLEMDATA::Traits::RangeType sigma_infty_grad_u_infty_chi;
sigma_infty_grad_u_infty *= (1 - param.get_chi(qp.position()));
sigma_infty_grad_u_infty *= (1.0 - param.get_chi(qp.position()));
RF factor = qp.weight() * geometry.integrationElement(qp.position());
for (std::size_t i = 0; i < lfsv.size(); i++)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment