Skip to content
Snippets Groups Projects
Commit bc5052be authored by Hörl, Maximilian's avatar Hörl, Maximilian
Browse files

implement penalty paramter for dg

parent 72f75c64
No related branches found
No related tags found
No related merge requests found
......@@ -26,6 +26,7 @@ public:
static constexpr auto edge = Dune::GeometryTypes::simplex(dim-1);
using Scalar = typename GridView::ctype;
using Coordinate = typename Problem::CoordinateType;
// using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,1,1>>;
// using Vector = Dune::BlockVector<Dune::FieldVector<Scalar,1>>;
......@@ -56,10 +57,10 @@ public:
// d = Vector(dof_);
}
void operator() (const Scalar mu)
void operator() (const Scalar mu0)
{
//assemble stiffness matrix A and load vector b
assembleSLE(mu);
assembleSLE(mu0);
//A->compress();
(*A).solve(d, b);
......@@ -143,7 +144,7 @@ protected:
}
//assemble stiffness matrix A and load vector b
void assembleSLE (const Scalar mu)
void assembleSLE (const Scalar mu0)
{
//quadrature rules
const QRuleElem& rule_q =
......@@ -166,6 +167,7 @@ protected:
{
const int elemIdx = mapper_.index(elem);
const auto& geo = elem.geometry();
const auto& center = geo.center();
//in the system of linear equations (SLE) Ad = b,
//the index elemIdxSLE refers to the basis function phi_elem,0
......@@ -233,6 +235,9 @@ protected:
const double intersctVol = intersctGeo.volume();
const auto& intersctCenter = intersctGeo.center();
//determine penalty parameter
const Scalar mu = getPenaltyParameter(mu0, intersct, center);
//create storage for multiply used integral values
//linearIntegrals[i] = int_intersct x[i] ds
Dune::FieldVector<Scalar, dim> linearIntegrals(0.0);
......@@ -572,6 +577,46 @@ private:
writeVTKoutput(dof_);
}
//returns the penalty parameter mu for a facet given as intersection of the
//grid where center is the center of intersection.inside()
const Scalar getPenaltyParameter (const Scalar mu0, const auto& intersection,
const auto& center) const
{
Scalar mu = getLargestEigenvalue(problem_.K(center))
/ getDiameter(intersection.inside());
if (intersection.neighbor())
{
const auto& neighbor = intersection.outside();
mu = std::max(mu,
getLargestEigenvalue(problem_.K(neighbor.geometry().center()))
/ getDiameter(neighbor));
}
return mu * mu0 * 2.0 * (1.0 + dim);
}
//returns the diameter of a grid element elem
const Scalar getDiameter (const auto& elem) const
{
Scalar diameter(0.0);
for (const auto& intersct : intersections(gridView_, elem))
{
diameter = std::max(diameter, intersct.geometry().volume());
}
return diameter;
}
//determines the largest absolute eigenvalue of a matrix
const Scalar getLargestEigenvalue (const auto& matrix) const
{
Coordinate eigenvalues;
Dune::FMatrixHelp::eigenValues(matrix, eigenvalues);
return eigenvalues.infinity_norm();
}
};
#endif
......@@ -9,6 +9,7 @@ class DGProblem
public:
static constexpr int dim = Coordinate::dimension;
using Matrix = Dune::FieldMatrix<Scalar, dim, dim>;
using CoordinateType = Coordinate;
//the exact solution at position pos
virtual Scalar exactSolution (const Coordinate& pos) const
......
......@@ -39,9 +39,7 @@ int main(int argc, char** argv)
//use YaspGrid if dim = 1 and MMesh otherwise
using Grid = std::conditional<dim == 1, Dune::YaspGrid<dim>,
Dune::MovingMesh<dim>>::type;
//using GridFactory = Dune::DGFGridFactory<Grid>;
using GridFactory = Dune::GmshGridFactory<Grid, /*implicit=*/false>;
using GridFactory = Dune::DGFGridFactory<Grid>;
using GridView = typename Grid::LeafGridView;
using Mapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
using Coordinate = Dune::FieldVector<double, dim>;
......@@ -53,11 +51,10 @@ int main(int argc, char** argv)
Dune::ParameterTreeParser::readINITree("parameterDG.ini", pt);
//assume constant penalty parameter
const double mu = std::stod(pt["mu"]);
const double mu0 = std::stod(pt["mu0"]);
//create a grid from .dgf file
//GridFactory gridFactory( "grids/cube" + std::to_string(dim) + "d.dgf" );
GridFactory gridFactory( "grids/mmdg2.msh" );
GridFactory gridFactory( "grids/cube" + std::to_string(dim) + "d.dgf" );
const Grid& grid = *gridFactory.grid();
const GridView& gridView = grid.leafGridView();
......@@ -89,7 +86,7 @@ int main(int argc, char** argv)
DG dg(gridView, mapper, *problem);
dg(mu);
dg(mu0);
//error and total run time
const double error = dg.computeL2error();
......
......@@ -8,7 +8,7 @@ from math import log
# writes a dgf file for a grid with N elements per directions and dimension dim
def writeDGF(N, dim):
# open dgf file
file = open(join("grids", "cube" + str(dim) + "d.dgf"), "rw")
file = open(join("grids", "cube" + str(dim) + "d.dgf"), "w")
file.write("DGF\n\nInterval\n")
......
......@@ -54,7 +54,7 @@ int main(int argc, char** argv)
Dune::ParameterTreeParser::readINITree("parameterMMDG.ini", pt);
//assume constant penalty parameter //TODO
const double mu = std::stod(pt["mu"]);
const double mu0 = std::stod(pt["mu0"]);
const double xi = std::stod(pt["xi"]); //coupling parameter
const double aperture = std::stod(pt["d"]);
......@@ -112,7 +112,7 @@ int main(int argc, char** argv)
const IMapper iMapper(iGridView, Dune::mcmgElementLayout());
MMDG mmdg(gridView, mapper, iGridView, iMapper, *problem);
mmdg(mu, xi);
mmdg(mu0, xi);
//error and total run time
const double error = mmdg.computeL2error();
......
mu = 1000
mu0 = 1000
problem = poisson
mu = 1000
mu0 = 1000
xi = 0.75
d = 1
problem = mmdg2
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment