Commit 3b2a72ac authored by Thomas Fetzer's avatar Thomas Fetzer
Browse files

[io]

adapted the interfacegridcreator to the common GridCreator type structure
and for the use in the common start.hh
deprecated the interfacemeshcreator.hh

[multidomain]
adapted multidomain to be used with the common start.hh
deprecated the old multidomainproblem constructor

thanks to and reviewed by bernd and gruenich



git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@15161 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 4bf4434e
......@@ -11,6 +11,10 @@ Differences Between DuMuX 2.7 and DuMuX 2.8
to the freeflow folder.
- Tests for coupling a turbulent free flow using zeroeq turbulence models
with flow in a porous medium have been added.
- The constructor of the multidomain problems has changed. They will now be called
directly from the start.hh, identical to the other problems. The new parameter
are the TimeManager and the HostGrid.
- The interfacegridcreator is now available as a new GridCreator type.
* IMMEDIATE INTERFACE CHANGES not allowing/requiring a deprecation period:
- The use and support for SGrid is dropped. SGrid is deprecated in Dune 2.4.
......@@ -22,10 +26,18 @@ Differences Between DuMuX 2.7 and DuMuX 2.8
the corresponding run.
* Deprecated CLASSES/FILES, to be removed after 2.8:
- The InterfaceMeshCreator has been moved to InterfaceGridCreator and adapted
to the structure of other grid creators, it can simply be used by specifying
the GridCreator TypeTag.
* Deprecated MEMBER FUNCTIONS, to be removed after 2.8:
- The method simulate(Scalar dtInitial, Scalar tEnd) from MultiDomainProblem,
it is unused and will be dropped.
- The GnuplotInterface functions are now called without giving a window number.
If plots should be plotted in different windows, different GnuplotInterface
objects are now required. This affects also all other plots in the "io" folder.
- The write() function in plotoverline2d.hh now has an append function, to be
able to decide whether the previously written file should be kept or not.
* Deprecated protected MEMBER VARIABLES, to be removed after 2.8: BEWARE: Older
compilers will not print any warning if a deprecated protected member variable
......
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
*
* \brief A grid creator for the coupled problems, with a refined interface
*/
#ifndef DUMUX_INTERFACEMESHCREATOR_HH
#define DUMUX_INTERFACEMESHCREATOR_HH
#include <dune/common/deprecated.hh>
#include <dune/grid/common/gridfactory.hh>
#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
#include <dune/grid/utility/structuredgridfactory.hh>
#include <dumux/common/basicproperties.hh>
namespace Dumux
{
namespace Properties
{
NEW_PROP_TAG(Scalar);
NEW_PROP_TAG(Grid);
}
/*!
* \brief A grid creator for the coupled problems, with a refined interface
*
* A grid creator, which can refine the grid towards the
* coupling interface and the top of the domain.
*/
template <class TypeTag>
class InterfaceGridCreator
{
public:
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef std::shared_ptr<Grid> GridPointer;
enum {dim = Grid::dimension};
/*!
* \brief Create the Grid
*/
static void makeGrid()
{
if (dim != 2)
{
DUNE_THROW(Dune::NotImplemented, "The InterfaceGridCreator is not implemented for 1D and 3D.");
}
Dune::array<unsigned int, dim> numCellsDummy = {{1,1}};
Dune::array<unsigned int, dim> numCells;
Dune::FieldVector<Scalar, dim> lowerLeft;
Dune::FieldVector<Scalar, dim> upperRight;
Dune::FieldVector<Scalar, dim> refinePoint(0);
Dune::FieldVector<Scalar, dim> gradingFactor(1);
// x-direction
numCells[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, unsigned int, Grid, NumberOfCellsX);
lowerLeft[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, LowerLeftX);
upperRight[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, UpperRightX);
// y-direction
numCells[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, unsigned int, Grid, NumberOfCellsY);
lowerLeft[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, LowerLeftY);
upperRight[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, UpperRightY);
refinePoint[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, InterfacePosY);
gradingFactor[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, GradingFactorY);
bool refineTop = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool, Grid, RefineTop);
typedef Dune::YaspGrid<dim> HelperGrid;
std::shared_ptr<HelperGrid> helperGrid = std::shared_ptr<HelperGrid> (
Dune::StructuredGridFactory<HelperGrid>::createCubeGrid(lowerLeft, upperRight, numCellsDummy));
typedef typename HelperGrid::LeafGridView HelperGridView;
typedef typename HelperGridView::template Codim<0>::Iterator HelperElementIterator;
typedef typename HelperGridView::Traits::template Codim<0>::Entity HelperElement;
typedef typename HelperElement::Geometry HelperGeometry;
HelperElementIterator helperElementIterator = helperGrid->leafGridView().template begin<0>();
const HelperElement& helperElement = *helperElementIterator;
const HelperGeometry& helperGeometry = helperElement.geometry();
Dune::FieldVector<Scalar,dim> refinePointLocal(helperGeometry.local(refinePoint));
std::cout << "refinePointGlobal = " << refinePoint
<< ", refinePointLocal = " << refinePointLocal
<< std::endl;
Dune::GridFactory<Grid> factory;
int nX = numCells[0];
int nY = numCells[1];
std::vector<std::vector<Scalar> > localPositions(dim);
for (int comp = 0; comp < dim; comp++)
{
Scalar lengthLeft = refinePointLocal[comp];
Scalar lengthRight = 1.0 - lengthLeft;
int nLeft, nRight;
Scalar hLeft, hRight;
if (lengthLeft < 1e-10)
{
nLeft = 0;
nRight = numCells[comp];
if (gradingFactor[comp] > 1.0)
hLeft = hRight = (1.0 - gradingFactor[comp])/(1.0 - pow(gradingFactor[comp], nRight));
else
hLeft = hRight = 1.0/numCells[comp];
}
else if (lengthLeft > 1.0 - 1e-10)
{
nLeft = numCells[comp];
nRight = 0;
if (gradingFactor[comp] > 1.0)
hLeft = hRight = (1.0 - gradingFactor[comp])/(1.0 - pow(gradingFactor[comp], nLeft));
else
hLeft = hRight = 1.0/numCells[comp];
}
else if (comp == dim - 1 && refineTop)
{
lengthLeft = refinePointLocal[comp];
lengthRight = (1 - refinePointLocal[comp])/2;
nLeft = nRight = numCells[comp]/3;
if (numCells[comp]%3 == 1)
nLeft += 1;
else if (numCells[comp]%3 == 2)
nRight += 1;
hLeft = lengthLeft*(1.0 - gradingFactor[comp])/(1.0 - pow(gradingFactor[comp], nLeft));
hRight = lengthRight*(1.0 - gradingFactor[comp])/(1.0 - pow(gradingFactor[comp], nRight));
}
else if (lengthLeft > 0.5)
{
Scalar nLeftDouble = std::ceil(-log((1.0 + sqrt(1.0 + 4.0 * pow(gradingFactor[comp], numCells[comp])
* lengthRight/lengthLeft))
/(2.0*pow(gradingFactor[comp], numCells[comp])))/log(gradingFactor[comp]));
nLeft = std::min((unsigned int)std::ceil(nLeftDouble), numCells[comp]);
nRight = numCells[comp] - nLeft;
if (gradingFactor[comp] > 1.0)
{
hLeft = lengthLeft*(1.0 - gradingFactor[comp])/(1.0 - pow(gradingFactor[comp], nLeft));
hRight = lengthRight*(1.0 - gradingFactor[comp])/(1.0 - pow(gradingFactor[comp], nRight));
}
else
hLeft = hRight = 1.0/numCells[comp];
}
else
{
Scalar nRightDouble = -log((1.0 + sqrt(1.0 + 4.0 * pow(gradingFactor[comp], numCells[comp])
* lengthLeft/lengthRight))
/(2.0*pow(gradingFactor[comp], numCells[comp])))/log(gradingFactor[comp]);
nRight = std::min((unsigned int)std::ceil(nRightDouble), numCells[comp]);
nLeft = numCells[comp] - nRight;
if (gradingFactor[comp] > 1.0)
{
hLeft = lengthLeft*(1.0 - gradingFactor[comp])/(1.0 - pow(gradingFactor[comp], nLeft));
hRight = lengthRight*(1.0 - gradingFactor[comp])/(1.0 - pow(gradingFactor[comp], nRight));
}
else
hLeft = hRight = 1.0/numCells[comp];
}
std::cout << "lengthLeft = " << lengthLeft
<< ", lengthRight = " << lengthRight
<< ", hLeft = " << hLeft
<< ", hRight = " << hRight
<< ", nLeft = " << nLeft
<< ", nRight = " << nRight
<< std::endl;
int numVertices = numCells[comp] + 1;
localPositions[comp].resize(numVertices);
localPositions[comp][0] = 0.0;
for (int i = 0; i < nLeft; i++)
{
Scalar hI = hLeft*pow(gradingFactor[comp], nLeft-1-i);
localPositions[comp][i+1] = localPositions[comp][i] + hI;
}
for (int i = 0; i < nRight; i++)
{
Scalar hI = hRight*pow(gradingFactor[comp], i);
localPositions[comp][nLeft+i+1] = localPositions[comp][nLeft+i] + hI;
}
if (comp == dim - 1 && refineTop)
for (int i = 0; i < nRight; i++)
{
Scalar hI = hRight*pow(gradingFactor[comp], nRight-1-i);
localPositions[comp][nLeft+nRight+i+1] = localPositions[comp][nLeft+nRight+i] + hI;
}
if (localPositions[comp][numVertices-1] != 1.0)
{
for (int i = 0; i < numVertices; i++)
localPositions[comp][i] /= localPositions[comp][numVertices-1];
}
}
Dune::FieldVector<Scalar,dim> local;
for (int j = 0; j < nY + 1; j++)
{
local[1] = localPositions[1][j];
for (int i = 0; i < nX + 1; i++)
{
local[0] = localPositions[0][i];
Dune::FieldVector<Scalar,dim> position(helperGeometry.global(local));
factory.insertVertex(position);
}
}
for (int j = 0; j < nY; j++)
{
for (int i = 0; i < nX; i++)
{
std::vector<unsigned int> vertices(4);
vertices[0] = j*(nX+1) + i;
vertices[1] = j*(nX+1) + i+1;
vertices[2] = (j+1)*(nX+1) + i;
vertices[3] = (j+1)*(nX+1) + i+1;
factory.insertElement(Dune::GeometryType(Dune::GeometryType::cube,dim), vertices);
}
}
gridPtr().reset(factory.createGrid());
}
DUNE_DEPRECATED_MSG("create() is deprecated use makeGrid() instead")
static Grid* create(const std::string& dgfName, const Dune::FieldVector<int, dim>& numCells,
const Scalar interfaceY, const Scalar grading, const bool refineTop = false)
{
typedef Dune::YaspGrid<dim> HelperGrid;
Dune::GridPtr<HelperGrid> helperGridPtr(dgfName);
HelperGrid& helperGrid = *helperGridPtr;
typedef typename HelperGrid::LeafGridView HelperGridView;
typedef typename HelperGridView::template Codim<0>::Iterator HelperElementIterator;
typedef typename HelperGridView::Traits::template Codim<0>::Entity HelperElement;
typedef typename HelperElement::Geometry HelperGeometry;
HelperElementIterator helperElementIterator = helperGrid.leafGridView().template begin<0>();
const HelperElement& helperElement = *helperElementIterator;
const HelperGeometry& helperGeometry = helperElement.geometry();
Dune::FieldVector<Scalar,dim> refinePoint(0);
refinePoint[dim-1] = interfaceY;
Dune::FieldVector<Scalar,dim> gradingFactor(1);
gradingFactor[dim-1] = grading;
Dune::FieldVector<Scalar,dim> refinePointLocal(helperGeometry.local(refinePoint));
std::cout << "rglobal = " << refinePoint
<< ", rlocal = " << refinePointLocal
<< std::endl;
Dune::GridFactory<Grid> factory;
int nX = numCells[0];
int nY = numCells[1];
std::vector<std::vector<Scalar> > localPositions(dim);
for (int comp = 0; comp < dim; comp++)
{
Scalar lengthLeft = refinePointLocal[comp];
Scalar lengthRight = 1.0 - lengthLeft;
int nLeft, nRight;
Scalar hLeft, hRight;
if (lengthLeft < 1e-10)
{
nLeft = 0;
nRight = numCells[comp];
if (gradingFactor[comp] > 1.0)
hLeft = hRight = (1.0 - gradingFactor[comp])/(1.0 - pow(gradingFactor[comp], nRight));
else
hLeft = hRight = 1.0/numCells[comp];
}
else if (lengthLeft > 1.0 - 1e-10)
{
nLeft = numCells[comp];
nRight = 0;
if (gradingFactor[comp] > 1.0)
hLeft = hRight = (1.0 - gradingFactor[comp])/(1.0 - pow(gradingFactor[comp], nLeft));
else
hLeft = hRight = 1.0/numCells[comp];
}
else if (comp == dim - 1 && refineTop)
{
lengthLeft = refinePointLocal[comp];
lengthRight = (1 - refinePointLocal[comp])/2;
nLeft = nRight = numCells[comp]/3;
if (numCells[comp]%3 == 1)
nLeft += 1;
else if (numCells[comp]%3 == 2)
nRight += 1;
hLeft = lengthLeft*(1.0 - gradingFactor[comp])/(1.0 - pow(gradingFactor[comp], nLeft));
hRight = lengthRight*(1.0 - gradingFactor[comp])/(1.0 - pow(gradingFactor[comp], nRight));
}
else if (lengthLeft > 0.5)
{
Scalar nLeftDouble = std::ceil(-log((1.0 + sqrt(1.0 + 4.0 * pow(gradingFactor[comp], numCells[comp])
* lengthRight/lengthLeft))
/(2.0*pow(gradingFactor[comp], numCells[comp])))/log(gradingFactor[comp]));
nLeft = std::min((int)std::ceil(nLeftDouble), numCells[comp]);
nRight = numCells[comp] - nLeft;
if (gradingFactor[comp] > 1.0)
{
hLeft = lengthLeft*(1.0 - gradingFactor[comp])/(1.0 - pow(gradingFactor[comp], nLeft));
hRight = lengthRight*(1.0 - gradingFactor[comp])/(1.0 - pow(gradingFactor[comp], nRight));
}
else
hLeft = hRight = 1.0/numCells[comp];
}
else
{
Scalar nRightDouble = -log((1.0 + sqrt(1.0 + 4.0 * pow(gradingFactor[comp], numCells[comp])
* lengthLeft/lengthRight))
/(2.0*pow(gradingFactor[comp], numCells[comp])))/log(gradingFactor[comp]);
nRight = std::min((int)std::ceil(nRightDouble), numCells[comp]);
nLeft = numCells[comp] - nRight;
if (gradingFactor[comp] > 1.0)
{
hLeft = lengthLeft*(1.0 - gradingFactor[comp])/(1.0 - pow(gradingFactor[comp], nLeft));
hRight = lengthRight*(1.0 - gradingFactor[comp])/(1.0 - pow(gradingFactor[comp], nRight));
}
else
hLeft = hRight = 1.0/numCells[comp];
}
std::cout << "lengthLeft = " << lengthLeft
<< ", lengthRight = " << lengthRight
<< ", hLeft = " << hLeft
<< ", hRight = " << hRight
<< ", nLeft = " << nLeft
<< ", nRight = " << nRight
<< std::endl;
int numVertices = numCells[comp] + 1;
localPositions[comp].resize(numVertices);
localPositions[comp][0] = 0.0;
for (int i = 0; i < nLeft; i++)
{
Scalar hI = hLeft*pow(gradingFactor[comp], nLeft-1-i);
localPositions[comp][i+1] = localPositions[comp][i] + hI;
}
for (int i = 0; i < nRight; i++)
{
Scalar hI = hRight*pow(gradingFactor[comp], i);
localPositions[comp][nLeft+i+1] = localPositions[comp][nLeft+i] + hI;
}
if (comp == dim - 1 && refineTop)
for (int i = 0; i < nRight; i++)
{
Scalar hI = hRight*pow(gradingFactor[comp], nRight-1-i);
localPositions[comp][nLeft+nRight+i+1] = localPositions[comp][nLeft+nRight+i] + hI;
}
if (localPositions[comp][numVertices-1] != 1.0)
{
for (int i = 0; i < numVertices; i++)
localPositions[comp][i] /= localPositions[comp][numVertices-1];
}
}
Dune::FieldVector<Scalar,dim> local;
for (int j = 0; j < nY + 1; j++)
{
local[1] = localPositions[1][j];
for (int i = 0; i < nX + 1; i++)
{
local[0] = localPositions[0][i];
Dune::FieldVector<Scalar,dim> position(helperGeometry.global(local));
factory.insertVertex(position);
}
}
for (int j = 0; j < nY; j++)
{
for (int i = 0; i < nX; i++)
{
std::vector<unsigned int> vertices(4);
vertices[0] = j*(nX+1) + i;
vertices[1] = j*(nX+1) + i+1;
vertices[2] = (j+1)*(nX+1) + i;
vertices[3] = (j+1)*(nX+1) + i+1;
factory.insertElement(Dune::GeometryType(Dune::GeometryType::cube,dim), vertices);
}
}
return factory.createGrid();
}
/*!
* \brief Returns a reference to the grid.
*/
static Grid &grid()
{
return *gridPtr();
}
/*!
* \brief Distributes the grid on all processes of a parallel
* computation.
*/
static void loadBalance()
{
gridPtr()->loadBalance();
}
/*!
* \brief Returns a reference to the shared pointer to the grid.
*/
static GridPointer &gridPtr()
{
static GridPointer interfaceGridCreator;
return interfaceGridCreator;
}
};
}
#endif
......@@ -28,6 +28,8 @@
#include <dune/grid/common/gridfactory.hh>
#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
#include <dune/common/deprecated.hh>
namespace Dumux
{
/*!
......@@ -38,7 +40,9 @@ namespace Dumux
*/
template <class Grid>
class InterfaceMeshCreator
class
DUNE_DEPRECATED_MSG("The InterfaceMeshCreator class is deprecated. Please use the the InterfaceGridCreator class instead.")
InterfaceMeshCreator
{
public:
enum {dim = Grid::dimension};
......
......@@ -56,6 +56,7 @@ private:
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
typedef typename GET_PROP_TYPE(TypeTag, GridCreator) GridCreator;
typedef typename GET_PROP_TYPE(TypeTag, SubDomain1TypeTag) SubDomain1TypeTag;
typedef typename GET_PROP_TYPE(TypeTag, SubDomain2TypeTag) SubDomain2TypeTag;
......@@ -69,8 +70,7 @@ private:
typedef typename GET_PROP_TYPE(SubDomain1TypeTag, GridView) SubDomainGridView1;
typedef typename GET_PROP_TYPE(SubDomain2TypeTag, GridView) SubDomainGridView2;
typedef typename GET_PROP_TYPE(TypeTag, MultiDomainGrid) MultiDomainGrid;
typedef typename GET_PROP_TYPE(TypeTag, MultiDomainGrid) MultiDomainGrid;
typedef typename MultiDomainGrid::LeafGridView MultiDomainGridView;
typedef typename MultiDomainGrid::Traits::template Codim<0>::Entity MultiDomainElement;
typedef typename MultiDomainGrid::SubDomainGrid SubDomainGrid;
......@@ -84,23 +84,32 @@ public:
*
* \param mdGrid The multi domain grid
* \param timeManager The TimeManager which is used by the simulation
*
*/
DUNE_DEPRECATED_MSG("This constructor is deprecated, please use the new constructor, which works with the common start facility. This has to be changed in the *.cc file.")
MultiDomainProblem(MultiDomainGrid &mdGrid,
TimeManager &timeManager)
: MultiDomainProblem(timeManager, GridCreator::grid().leafGridView())
{}
/*!
* \brief The problem for the coupling of Stokes and Darcy flow
*
* \param timeManager The time manager
* \param gridView The grid view
*/
template<class GridView>
MultiDomainProblem(TimeManager &timeManager,
GridView gridView)
: timeManager_(timeManager)
, newtonMethod_(asImp_())
, newtonCtl_(asImp_())
, mdGrid_(mdGrid)
, mdGridView_(mdGrid.leafGridView())
, mdVertexMapper_(mdGrid_.leafGridView())
, sdID1_(0)
, sdID2_(1)
, sdGrid1_(mdGrid.subDomain(sdID1_))
, sdGrid2_(mdGrid.subDomain(sdID2_))
, sdProblem1_(timeManager, sdGrid1_.leafGridView())
, sdProblem2_(timeManager, sdGrid2_.leafGridView())
{}
{
mdGrid_ = std::make_shared<MultiDomainGrid> (GridCreator::grid());
mdGridView_ = std::make_shared<MultiDomainGridView> (mdGrid_->leafGridView());
mdVertexMapper_ = std::make_shared<VertexMapper> (mdGrid_->leafGridView());
sdProblem1_ = std::make_shared<SubDomainProblem1> (timeManager, mdGrid_->subDomain(sdID1()).leafGridView());
sdProblem2_ = std::make_shared<SubDomainProblem2> (timeManager, mdGrid_->subDomain(sdID2()).leafGridView());
}
//! \copydoc Dumux::ImplicitProblem::init()
void init()
......@@ -112,7 +121,7 @@ public:
// set the initial condition of the model
model().init(asImp_());
// intialize lagrange multipliers
// initialize Lagrange multipliers
asImp_().initMortarElements();
}
......@@ -130,7 +139,6 @@ public:
<