diff --git a/dune/mmdg/problems/mmdgproblem4.hh b/dune/mmdg/problems/mmdgproblem4.hh index cf09ef6b727e7b946854d364c29605e582e9b6ba..347d4096782730bf0c66a586636b4a8878b6d19f 100644 --- a/dune/mmdg/problems/mmdgproblem4.hh +++ b/dune/mmdg/problems/mmdgproblem4.hh @@ -10,6 +10,7 @@ class MMDGProblem4 : public CoupledDGProblem<Coordinate, Scalar> { public: static constexpr int dim = Coordinate::dimension; + static constexpr Scalar jump = 1.0; using Base = CoupledDGProblem<Coordinate, Scalar>; using Matrix = typename Base::Matrix; @@ -23,7 +24,7 @@ public: if (pos[0] > 0.5) { - return solution += d_; + return solution += jump; } return solution; @@ -32,7 +33,7 @@ public: //the exact solution on the interface at position pos Scalar exactInterfaceSolution (const Coordinate& pos) const { - return 0.25 - pos[1] * pos[1] + 0.5 * d_; + return 0.25 - pos[1] * pos[1] + 0.5 * jump; } //indicates whether an exact solution is implemented for the problem @@ -44,19 +45,34 @@ public: //interface source term at position pos Scalar qInterface (const Coordinate& pos) const { - return 2.0 * d_; + // return 2.0 * d_ * d_; + const Scalar de = aperture(pos); + const Scalar dPrime = sqrt(gradAperture(pos) * gradAperture(pos)); + const Scalar dPrimePrime = 0.; + + return (d_ * d_ * ( 10. * pos[1] * 0.01 + 9. * pos[1] * pos[1] + 2. * 0.01 * 0.01 - 0.5 * (jump + 0.5))) * de; } //aperture d of the fracture at position pos Scalar aperture (const Coordinate& pos) const { - return d_; + return d_ * (pos[1] + 0.01);; + } + + //tangential gradient of the aperture d of the fracture at position pos + Coordinate gradAperture (const Coordinate& pos) const + { + Coordinate gradD(0.0); + gradD[1] = d_; + return gradD; + + return Coordinate(0.0); } //permeability of the fracture in normal direction at position pos Scalar Kperp (const Coordinate& pos) const { - return Scalar(1.0) / d_; + return Scalar(1.0) / jump; } //tangential permeability tensor of the interface at position pos @@ -66,7 +82,7 @@ public: for (int i = 0; i < dim; i++) { - permeability[i][i] = 1.0 / d_; + permeability[i][i] = 1.0 * aperture(pos);//d_; } return permeability; @@ -76,7 +92,7 @@ public: //over x * boundary(x) int quadratureOrderBoundary () const { - return 3; + return 10;//3; } //returns the recommended quadrature order to compute an integral @@ -85,7 +101,7 @@ public: //and d^2 * interfaceBoundary(x) * Kpar int quadratureOrderInterfaceBoundary () const { - return 3; + return 10;//3; } private: