diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh index dda8217a47f3c2a6ccf67adac688e63cbafd125b..24a4c4b48c4501ded175ee59a29103037833869c 100644 --- a/dune/mmdg/mmdg.hh +++ b/dune/mmdg/mmdg.hh @@ -28,9 +28,9 @@ public: (1 + dim) * gridView.size(0) + dim * iGridView.size(0)), iGridView_(iGridView), iMapper_(iMapper) {} - const void operator() (const Scalar mu) + const void operator() (const Scalar mu, const Scalar xi) { - assembleSLE(mu); + assembleSLE(mu, xi); Base::A->compress(); } @@ -39,11 +39,28 @@ public: private: //assemble stiffness matrix A and load vector b - void assembleSLE (const Scalar mu) + void assembleSLE (const Scalar mu, const Scalar xi) { + //modell parameters //NOTE: what should be the capture of alpha? + const auto alpha = [this](const Coordinate& pos) -> Scalar + { + return Base::problem_.Kperp(pos) / Base::problem_.aperture(pos); + }; + + const auto beta = [&alpha, &xi](const Coordinate& pos) -> Scalar + { + return 4 * alpha(pos) / (2 * xi - 1); + }; + //quadrature rules const QRuleInterface& rule_qInterface = QRulesInterface::rule(iSimplex, Base::problem_.quadratureOrder_qInterface()); + const QRuleInterface& rule_Kpar = QRulesInterface::rule(iSimplex, + Base::problem_.quadratureOrder_Kparallel()); + const QRuleInterface& rule_Kperp = QRulesInterface::rule(iSimplex, + Base::problem_.quadratureOrder_Kperp()); + const QRuleInterface& rule_KperpPlus1 = QRulesInterface::rule(iSimplex, + Base::problem_.quadratureOrder_Kperp() + 1); //we use the bulk basis // phi_elem,0 (x) = indicator(elem); @@ -75,10 +92,11 @@ private: { const int iElemIdx = iMapper_.index(iElem); const auto& iGeo = iElem.geometry(); + const auto& intersct = Base::grid_.asIntersection(iElem); //get basis of the interface coordinate system for evaluation of //basis function phi_iElem,i with i=1,...,dim-1 - const std::array<Coordinate, dim-1> iFrame = interfaceFrame( + const InterfaceBasis iFrame = interfaceFrame( Base::grid_.asIntersection(iElem).centerUnitOuterNormal()); //in the total system of linear equations (SLE) Ad = b, @@ -111,11 +129,80 @@ private: //TODO: fill stiffness matrix A with interface entries - //fill stiffness matrix A with coupling entries - for (const auto& elem : elements(Base::gridView_)) + //evaluate + // int_iElem d * (K_par*grad_par(phi_iElem,i)) * grad_par(phi_iElem,j) ds + // = int_iElem d * (K_par * tau_i) * tau_j ds + //using a quadrature rule for i,j = 1,...,dim-1 + for (const auto& ip : rule_Kpar) + { + const auto& qpLocal = ip.position(); + const auto& qpGlobal = iGeo.global(qpLocal); + + const Scalar weight(ip.weight()); + const Scalar integrationElem(iGeo.integrationElement(qpLocal)); + + for (int i = 0; i < dim-1 ; i++) + { + Coordinate K_parDotTau_i; + Base::problem_.Kparallel(qpGlobal).mv(iFrame[i], K_parDotTau_i); + + for (int j = 0; j <= i; j++) //TODO: symmetrize + { + Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) += + weight * integrationElem * Base::problem_.aperture(qpGlobal) + * ( K_parDotTau_i * iFrame[j] ); + } + } + } + + //TODO: check interface bilinear form + + //evaluate the coupling term + // int_iElem beta * phi_iElem,i * phi_iElem,j ds + //using a quadrature rule for i = 0 and j = 0,...,dim-1 //TODO or vice versa (symmetrize) + for (const auto& ip : rule_Kperp) + { + const auto& qpLocal = ip.position(); + const auto& qpGlobal = iGeo.global(qpLocal); + + const Scalar weight(ip.weight()); + const Scalar integrationElem(iGeo.integrationElement(qpLocal)); + + Base::A->entry(iElemIdxSLE, iElemIdxSLE) += + weight * integrationElem * beta(qpGloabal); + + for (int i = 0; i < dim - 1; i++) + { + Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) += + weight * integrationElem * beta(qpGlobal) * (iFrame[i] * qpGlobal) + * (iFrame[j] * qpGlobal); + } + } + + //evaluate the coupling term + // int_iElem beta * phi_iElem,i * phi_iElem,j ds + //using a quadrature rule for i,j = 1,...,dim-1 //TODO symmetrize + for (const auto& ip : rule_KperpPlus1) { - //TODO + const auto& qpLocal = ip.position(); + const auto& qpGlobal = iGeo.global(qpLocal); + + const Scalar weight(ip.weight()); + const Scalar integrationElem(iGeo.integrationElement(qpLocal)); + + for (int i = 0; i < dim - 1; i++) + { + for (int j = 0; j <= i; j++) + { + Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE) += + weight * integrationElem * beta(qpGlobal) * (iFrame[i] * qpGlobal) ; + } + } } + + const auto elem1 = intersct.inside(); + const auto elem2 = intersct.ouside(); + //TODO coupling entries } } diff --git a/dune/mmdg/problems/coupleddgproblem.hh b/dune/mmdg/problems/coupleddgproblem.hh index e1abfffe97541843449c9014154875354d757930..ad26400cc62d663d9948d01156fa6f963c736c30 100644 --- a/dune/mmdg/problems/coupleddgproblem.hh +++ b/dune/mmdg/problems/coupleddgproblem.hh @@ -20,12 +20,18 @@ class CoupledDGProblem : public DGProblem<Vector, Scalar> return Scalar(0.0); } - //aperture of the fracture at position pos + //aperture d of the fracture at position pos virtual Scalar aperture (const Vector& pos) const { return Scalar(1.0); } + //(tangential) gradient of the aperture at position pos + virtual Vector gradAperture (const Vector& pos) const + { + return Vector(0.0); + } + //tangential permeability tensor of the interface at position pos virtual Matrix Kparallel (const Vector& pos) const { @@ -53,6 +59,20 @@ class CoupledDGProblem : public DGProblem<Vector, Scalar> { return 1; } + + //returns the recommended quadrature order to compute an integral + //over d(x) * Kparallel(x) + virtual int quadratureOrder_Kparallel () const + { + return 1; + } + + //returns the recommended quadrature order to compute an integral + //over x * Kperp(x) / d(x) + virtual int quadratureOrder_Kperp () const + { + return 1; + } }; #endif diff --git a/src/mmdg.cc b/src/mmdg.cc index ed37ba1263c0c2d16bb5327c4fc9c7dde4ff481c..32140f3f014c8f9c21fb44d8be56ddeed52d7955 100644 --- a/src/mmdg.cc +++ b/src/mmdg.cc @@ -48,8 +48,9 @@ int main(int argc, char** argv) Dune::ParameterTree pt; Dune::ParameterTreeParser::readINITree("parameterMMDG.ini", pt); - //assume constant penalty parameter + //assume constant penalty parameter //TODO const double mu = std::stod(pt["mu"]); + const double xi = std::stod(pt["xi"]); //coupling parameter //create a grid from .dgf file GridFactory gridFactory( "grids/cube" + std::to_string(dim) + "d.dgf" ); @@ -77,7 +78,7 @@ int main(int argc, char** argv) } MMDG mmdg(grid, gridView, mapper, iGridView, iMapper, *problem); - mmdg(mu); + mmdg(mu, xi); //print total run time const double timeTotal = ClockHelper::elapsedRuntime(timeInitial); diff --git a/src/parameterMMDG.ini b/src/parameterMMDG.ini index 165571303cbed4b80a298d46af9b7007f682ed5d..c787684f0ec215e6247dd051e04207ae13516ea5 100644 --- a/src/parameterMMDG.ini +++ b/src/parameterMMDG.ini @@ -1,2 +1,3 @@ mu = 1000 +xi = 0.75 problem = poisson