Commit c443dcd2 authored by Timo Koch's avatar Timo Koch
Browse files

[gridcreator] Support gmsh physical entities. Add a test for gmsh

functionality.

Add a test for the gmsh interface of the gridcreator. The example uses
parametrized boundary segments read from the gmsh file and
boundaryMarkers defined in the gmsh file. The output files are fuzzy
comopared.

Reviewed by bernd


git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@15348 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 5bcb2486
......@@ -134,6 +134,38 @@ public:
DUNE_THROW(Dune::InvalidStateException, "The parameters method is only available if the grid was constructed with a DGF file!");
}
/*!
* \brief Call the parameters function of the DGF grid pointer if available
*/
static const int getBoundaryDomainMarker(int boundarySegmentIndex)
{
if(enableGmshDomainMarkers_)
return boundaryMarkers_[boundarySegmentIndex];
else
DUNE_THROW(Dune::InvalidStateException, "The getBoundaryDomainMarker method is only available if DomainMarkers for Gmsh were enabled!"
<< " If your Gmsh file contains domain markers / physical entities,"
<< " enable them by setting " << GET_PROP_VALUE(TypeTag, GridParameterGroup)
<< ".DomainMarkers = 1 in the input file.");
}
/*!
* \brief Call the parameters function of the DGF grid pointer if available
*/
static const int getElementDomainMarker(int elementIdx)
{
if(enableGmshDomainMarkers_)
{
if(elementIdx >= grid().levelGridView(0).size(0))
DUNE_THROW(Dune::RangeError, "Requested element index is bigger than the number of level 0 elements!");
return elementMarkers_[elementIdx];
}
else
DUNE_THROW(Dune::InvalidStateException, "The getElementDomainMarker method is only available if DomainMarkers for Gmsh were enabled!"
<< " If your Gmsh file contains domain markers / physical entities,"
<< " enable them by setting " << GET_PROP_VALUE(TypeTag, GridParameterGroup)
<< ".DomainMarkers = 1 in the input file.");
}
/*!
* \brief Call loadBalance() function of the grid.
*/
......@@ -162,7 +194,7 @@ protected:
}
/*!
* \brief Returns a reference to the dgf grid pointer (Dune::GridPtr<Grid>).
* \brief Returns a reference to the DGF grid pointer (Dune::GridPtr<Grid>).
*/
static Dune::GridPtr<Grid> &dgfGridPtr()
{
......@@ -172,7 +204,7 @@ protected:
return dgfGridPtr_;
}
else
DUNE_THROW(Dune::InvalidStateException, "The dgf grid pointer is only available if the grid was constructed with a DGF file!");
DUNE_THROW(Dune::InvalidStateException, "The DGF grid pointer is only available if the grid was constructed with a DGF file!");
}
/*!
......@@ -212,14 +244,26 @@ protected:
{
// get some optional parameters
bool verbose = false;
try { verbose = GET_RUNTIME_PARAM_FROM_GROUP_CSTRING(TypeTag, bool, GET_PROP_VALUE(TypeTag, GridParameterGroup).c_str(), Verbose);}
try { verbose = GET_RUNTIME_PARAM_FROM_GROUP_CSTRING(TypeTag, bool, GET_PROP_VALUE(TypeTag, GridParameterGroup).c_str(), Verbosity);}
catch (Dumux::ParameterException &e) { }
bool boundarySegments = false;
try { boundarySegments = GET_RUNTIME_PARAM_FROM_GROUP_CSTRING(TypeTag, bool, GET_PROP_VALUE(TypeTag, GridParameterGroup).c_str(), BoundarySegments);}
catch (Dumux::ParameterException &e) { }
gridPtr() = std::shared_ptr<Grid>(Dune::GmshReader<Grid>::read(fileName, verbose, boundarySegments));
bool domainMarkers = false;
try { domainMarkers = GET_RUNTIME_PARAM_FROM_GROUP_CSTRING(TypeTag, bool, GET_PROP_VALUE(TypeTag, GridParameterGroup).c_str(), DomainMarkers);}
catch (Dumux::ParameterException &e) { }
if(domainMarkers)
{
enableGmshDomainMarkers_ = true;
gridPtr() = std::shared_ptr<Grid>(Dune::GmshReader<Grid>::read(fileName, boundaryMarkers_, elementMarkers_, verbose, boundarySegments));
}
else
{
gridPtr() = std::shared_ptr<Grid>(Dune::GmshReader<Grid>::read(fileName, verbose, boundarySegments));
}
}
}
......@@ -290,11 +334,32 @@ protected:
* It is always enabled if a DGF grid file was used to create the grid.
*/
static bool enableDgfGridPointer_;
/*!
* \brief A state variable if domain markers have been read from a gmsh file.
*/
static bool enableGmshDomainMarkers_;
/*!
* \brief Element and domain markers obtained from Gmsh physical entities
* They map from element indices / boundary ids to the physical entity mark
*/
static std::vector<int> elementMarkers_;
static std::vector<int> boundaryMarkers_;
};
template <class TypeTag, class Grid>
bool GridCreatorBase<TypeTag, Grid>::enableDgfGridPointer_ = false;
template <class TypeTag, class Grid>
bool GridCreatorBase<TypeTag, Grid>::enableGmshDomainMarkers_ = false;
template <class TypeTag, class Grid>
std::vector<int> GridCreatorBase<TypeTag, Grid>::elementMarkers_;
template <class TypeTag, class Grid>
std::vector<int> GridCreatorBase<TypeTag, Grid>::boundaryMarkers_;
/*!
* \brief Provides the grid creator implementation for all supported grid managers that constructs a grid
* from information in the input file. This class is specialised below for all
......@@ -674,7 +739,7 @@ private:
* - Cells : number of elements in a structured grid
* - CellType : "Cube" or "Simplex" to be used for structured grids
* - Refinement : the number of global refines to perform
* - Verbose : whether the grid construction should output to standard out
* - Verbosity : whether the grid construction should output to standard out
* - HeapSize: The heapsize used to allocate memory
* - BoundarySegments : whether to insert boundary segments into the grid
*
......@@ -708,7 +773,7 @@ public:
preProcessing_();
// Check for cell type
std::string cellType = "Cube";
try { cellType = GET_RUNTIME_PARAM_FROM_GROUP_CSTRING(TypeTag, int, GET_PROP_VALUE(TypeTag, GridParameterGroup).c_str(), CellType);
try { cellType = GET_RUNTIME_PARAM_FROM_GROUP_CSTRING(TypeTag, std::string, GET_PROP_VALUE(TypeTag, GridParameterGroup).c_str(), CellType);
if (cellType != "Cube" && cellType != "Simplex")
DUNE_THROW(Dune::IOError, "UGGrid only supports 'Cube' or 'Simplex' as closure type. Not '"<< cellType<<"'!");
}
......@@ -800,7 +865,7 @@ private:
* - UpperRight : upperright corner of a structured grid
* - Cells : number of elements in a structured grid
* - Refinement : the number of global refines to perform
* - Verbose : whether the grid construction should output to standard out
* - Verbosity : whether the grid construction should output to standard out
* - BoundarySegments : whether to insert boundary segments into the grid
*
*/
......@@ -891,7 +956,7 @@ public:
* The following keys are recognized:
* - File : A DGF or gmsh file to load from, type detection by file extension
* - Verbose : whether the grid construction should output to standard out
* - Verbosity : whether the grid construction should output to standard out
* - LowerLeft : lowerleft corner of a structured grid
* - UpperRight : upperright corner of a structured grid
* - Cells : number of elements in a structured grid
......@@ -943,7 +1008,7 @@ public:
* The following keys are recognized:
* - File : A DGF or gmsh file to load from, type detection by file extension
* - Verbose : whether the grid construction should output to standard out
* - Verbosity : whether the grid construction should output to standard out
* - LowerLeft : lowerleft corner of a structured grid
* - UpperRight : upperright corner of a structured grid
* - Cells : number of elements in a structured grid
......@@ -1028,7 +1093,7 @@ public:
* The following keys are recognized:
* - File : A DGF or gmsh file to load from, type detection by file extension
* - Verbose : whether the grid construction should output to standard out
* - Verbosity : whether the grid construction should output to standard out
*
*/
template<class TypeTag, int dim, int dimworld>
......
add_subdirectory("gnuplotinterface")
\ No newline at end of file
add_subdirectory("gnuplotinterface")
add_subdirectory("gridcreator")
\ No newline at end of file
add_dumux_test(test_gridcreator_gmsh test_gridcreator_gmsh test_gridcreator_gmsh.cc
python ${CMAKE_SOURCE_DIR}/bin/runtest.py
--script fuzzy
--command "${CMAKE_CURRENT_BINARY_DIR}/test_gridcreator_gmsh"
--files ${CMAKE_SOURCE_DIR}/test/references/bifurcation-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/bifurcation-00000.vtu
${CMAKE_SOURCE_DIR}/test/references/bifurcation-reference-refined.vtu
${CMAKE_CURRENT_BINARY_DIR}/bifurcation-00001.vtu)
size = 0.7;
length = 10.0;
radius = 1.0;
a = 30/360*2*Pi;
// quadratic elements for boundary parametrization
Mesh.ElementOrder = 2;
Point(1) = {0.0, 0.0, 0.0, size};
Point(17) = {0.0, 0.0, radius, size};
Point(18) = {0.0, 0.0, -radius, size};
Point(2) = {0.0, length, 0.0, size};
Point(3) = {radius, length, 0.0, size};
Point(4) = {-radius, length, 0.0, size};
Point(5) = {0.0, length, radius, size};
Point(6) = {0.0, length, -radius, size};
Point(7) = {length*Cos(a), -length*Sin(a), 0.0, size};
Point(8) = {length*Cos(a) + radius*Sin(a), -length*Sin(a) + radius*Cos(a), 0.0, size};
Point(9) = {length*Cos(a) - radius*Sin(a), -length*Sin(a) - radius*Cos(a), 0.0, size};
Point(10) = {length*Cos(a), -length*Sin(a), radius, size};
Point(11) = {length*Cos(a), -length*Sin(a), -radius, size};
Point(12) = {-length*Cos(a), -length*Sin(a), 0.0, size};
Point(13) = {-length*Cos(a) - radius*Sin(a), -length*Sin(a) + radius*Cos(a), 0.0, size};
Point(14) = {-length*Cos(a) + radius*Sin(a), -length*Sin(a) - radius*Cos(a), 0.0, size};
Point(15) = {-length*Cos(a), -length*Sin(a), radius, size};
Point(16) = {-length*Cos(a), -length*Sin(a), -radius, size};
Point(19) = {radius, 2.0*radius*Cos(a)/3.0, 0.0, size};
Point(20) = {-radius, 2.0*radius*Cos(a)/3.0, 0.0, size};
Point(21) = {0.0, -4.0*radius*Cos(a)/3.0, 0.0, size};
//optional circle segments
//Point(22) = {0, 2.0*radius*Cos(a)/3.0, 0.0, size};
//Point(23) = {0, 2.0*radius*Cos(a)/3.0, radius, size};
//Point(24) = {0, 2.0*radius*Cos(a)/3.0, -radius, size};
Circle(1) = {15, 12, 14};
Circle(2) = {14, 12, 16};
Circle(3) = {16, 12, 13};
Circle(4) = {13, 12, 15};
Circle(5) = {4, 2, 6};
Circle(6) = {4, 2, 5};
Circle(7) = {5, 2, 3};
Circle(8) = {6, 2, 3};
Circle(9) = {8, 7, 10};
Circle(10) = {10, 7, 9};
Circle(11) = {9, 7, 11};
Circle(12) = {11, 7, 8};
Line(13) = {15, 17};
Line(15) = {17, 10};
Line(17) = {17, 5};
Line(18) = {6, 18};
Line(19) = {18, 11};
Line(20) = {16, 18};
Line(21) = {13, 20};
Line(22) = {4, 20};
Line(23) = {3, 19};
Line(24) = {19, 8};
Line(25) = {21, 9};
Line(26) = {21, 14};
Ellipse(27) = {17, 1, 18, 21};
Ellipse(28) = {18, 1, 17, 21};
Ellipse(29) = {17, 1, 18, 19};
Ellipse(30) = {18, 1, 17, 19};
Ellipse(31) = {17, 1, 18, 20};
Ellipse(32) = {18, 1, 17, 20};
//Circle(33) = {24, 22, 19};
//Circle(34) = {19, 22, 23};
//Circle(35) = {23, 22, 20};
//Circle(36) = {20, 22, 24};
Line Loop(33) = {5, 8, -7, -6};
Plane Surface(34) = {33};
Line Loop(35) = {11, 12, 9, 10};
Plane Surface(36) = {35};
Line Loop(37) = {3, 4, 1, 2};
Plane Surface(38) = {37};
Line Loop(39) = {18, 32, -22, 5};
Ruled Surface(40) = {39};
Line Loop(41) = {32, -21, -3, 20};
Ruled Surface(42) = {41};
Line Loop(43) = {18, 30, -23, -8};
Ruled Surface(44) = {43};
Line Loop(45) = {22, -31, 17, -6};
Ruled Surface(46) = {45};
Line Loop(47) = {7, 23, -29, 17};
Ruled Surface(48) = {47};
Line Loop(49) = {21, -31, -13, -4};
Ruled Surface(50) = {49};
Line Loop(51) = {13, 27, 26, -1};
Ruled Surface(52) = {51};
Line Loop(53) = {20, 28, 26, 2};
Ruled Surface(54) = {53};
Line Loop(55) = {25, -10, -15, 27};
Ruled Surface(56) = {55};
Line Loop(57) = {25, 11, -19, 28};
Ruled Surface(58) = {57};
Line Loop(59) = {15, -9, -24, -29};
Ruled Surface(60) = {59};
Line Loop(61) = {19, 12, -24, -30};
Ruled Surface(62) = {61};
Line(64) = {17, 1};
Line(65) = {1, 18};
Line Loop(66) = {29, -30, -65, -64};
Plane Surface(67) = {66};
Line Loop(68) = {28, -27, 64, 65};
Plane Surface(69) = {68};
Line Loop(70) = {31, -32, -65, -64};
Plane Surface(71) = {70};
Surface Loop(72) = {34, 40, 44, 48, 46, 67, 71};
Volume(73) = {72};
Surface Loop(74) = {42, 50, 52, 54, 38, 71, 69};
Volume(75) = {74};
Surface Loop(76) = {58, 56, 36, 62, 60, 69, 67};
Volume(77) = {76};
Physical Surface(1) = {34};
Physical Surface(2) = {36};
Physical Surface(3) = {38};
Physical Surface(4) = {40, 42, 44, 46, 48, 50, 52, 54, 56, 58, 60, 62};
Physical Volume(1) = {73};
Physical Volume(2) = {75};
Physical Volume(3) = {77};
This diff is collapsed.
/*****************************************************************************
* 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 Test for gmsh interface of the grid creator
*/
#include "config.h"
#include <iostream>
#include <dune/common/parametertreeparser.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/grid/io/file/vtk.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dumux/io/gridcreator.hh>
#include <dumux/common/basicproperties.hh>
namespace Dumux {
template<class TypeTag>
class GridCreatorGmshTest;
namespace Properties
{
NEW_TYPE_TAG(GridCreatorGmshTest, INHERITS_FROM(NumericModel));
SET_TYPE_PROP(GridCreatorGmshTest, Grid, Dune::UGGrid<3>);
// Change the default "Grid" to customized "BifurcationGrid", merely for demonstration purposes.
SET_STRING_PROP(GridCreatorGmshTest, GridParameterGroup, "BifurcationGrid");
}
template<class TypeTag>
class GridCreatorGmshTest
{
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
static const int dim = Grid::dimension;
typedef typename Dumux::GridCreator<TypeTag> GridCreator;
typedef typename Dune::ReferenceElements<Scalar, dim> ReferenceElements;
typedef typename Dune::ReferenceElement<Scalar, dim> ReferenceElement;
typedef typename Dune::LeafMultipleCodimMultipleGeomTypeMapper<Grid, Dune::MCMGVertexLayout> VertexMapper;
public:
static void getBoundaryDomainMarkers(std::vector<int>& boundaryMarker)
{
const auto& gridView = GridCreator::grid().leafGridView();
VertexMapper vertexMapper(GridCreator::grid());
boundaryMarker.clear();
boundaryMarker.resize(gridView.size(dim));
for(auto eIt = gridView.template begin<0>(); eIt != gridView.template end<0>(); ++eIt)
{
for(auto isIt = gridView.ibegin(*eIt); isIt != gridView.iend(*eIt); ++isIt)
{
if(!isIt->boundary())
continue;
const ReferenceElement &refElement = ReferenceElements::general(eIt->geometry().type());
// loop over vertices of the intersection facet
for(int vIdx = 0; vIdx < refElement.size(isIt->indexInInside(), 1, dim); vIdx++)
{
// get local vertex index with respect to the element
int vIdxLocal = refElement.subEntity(isIt->indexInInside(), 1, vIdx, dim);
auto vertex = eIt->template subEntity<dim>(vIdxLocal);
// make sure we always take the lowest non-zero marker (problem dependent!)
if (boundaryMarker[vertexMapper.index(vertex)] == 0)
boundaryMarker[vertexMapper.index(vertex)] = GridCreator::getBoundaryDomainMarker(isIt->boundarySegmentIndex());
else
{
if (boundaryMarker[vertexMapper.index(vertex)] > GridCreator::getBoundaryDomainMarker(isIt->boundarySegmentIndex()))
boundaryMarker[vertexMapper.index(vertex)] = GridCreator::getBoundaryDomainMarker(isIt->boundarySegmentIndex());
}
}
}
}
}
};
}
int main(int argc, char** argv)
{
#if HAVE_UG
try {
// Some typedefs
typedef typename TTAG(GridCreatorGmshTest) TypeTag;
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename Dumux::GridCreator<TypeTag> GridCreator;
// Read the parameters from the input file
typedef typename GET_PROP(TypeTag, ParameterTree) ParameterTree;
Dune::ParameterTreeParser::readINITree("test_gridcreator_gmsh.input", ParameterTree::tree());
// Make the grid
GridCreator::makeGrid();
// Read the boundary markers and convert them to vertex flags (e.g. for use in a box method)
// Write a map from vertex position to boundaryMarker
std::vector<int> boundaryMarker;
Dumux::GridCreatorGmshTest<TypeTag>::getBoundaryDomainMarkers(boundaryMarker);
// construct a vtk output writer and attach the boundaryMakers
Dune::VTKSequenceWriter<Grid::LeafGridView> vtkWriter(GridCreator::grid().leafGridView(), "bifurcation", ".", "");
vtkWriter.addVertexData(boundaryMarker, "boundaryMarker");
vtkWriter.write(0);
// refine grid once. Due to parametrized boundaries this will result in a grid closer to the orginal geometry.
GridCreator::grid().globalRefine(1);
Dumux::GridCreatorGmshTest<TypeTag>::getBoundaryDomainMarkers(boundaryMarker);
vtkWriter.write(1);
}
catch (Dumux::ParameterException &e) {
typedef typename TTAG(GridCreatorGmshTest) TypeTag;
Dumux::Parameters::print<TypeTag>();
std::cerr << e << ". Abort!\n";
return 1;
}
catch (Dune::Exception &e) {
std::cerr << "Dune reported error: " << e << std::endl;
return 3;
}
catch (...) {
std::cerr << "Unknown exception thrown!\n";
return 4;
}
#else
#warning "You need to have UGGrid installed to run this test."
std::cerr << "You need to have UGGrid installed to run this test\n";
return 77;
#endif
}
[BifurcationGrid]
File = ./grids/bifurcation.msh
Verbosity = 1 # verbose grid file parsing
BoundarySegments = 1 # read parametrized boundary segments
DomainMarkers = 1 # read domain markers (gmsh physical entities)
This diff is collapsed.
This diff is collapsed.
Supports Markdown
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