Commit 45c1d4a7 authored by Timo Koch's avatar Timo Koch
Browse files

Merge branch 'fix/bugs-cakegridmanager' into 'master'

[cakegridmanager] Fix bugs, allow zero well radius

See merge request !1785
parents 7593ed2d 3945f949
......@@ -18,6 +18,7 @@ from six.moves import zip
import os
import re
import glob
import functools
# fuzzy compare VTK tree from VTK strings
def compare_vtk(vtk1, vtk2, absolute=1.5e-7, relative=1e-2, zeroValueThreshold={}, verbose=True):
......@@ -75,8 +76,10 @@ def compare_vtk(vtk1, vtk2, absolute=1.5e-7, relative=1e-2, zeroValueThreshold={
# do the fuzzy compare
if is_fuzzy_equal_node(sortedroot1, sortedroot2, absolute, relative, zeroValueThreshold, verbose, convertedFromParallelVtu):
print("Fuzzy comparison done (equal)")
return 0
else:
print("Fuzzy comparison done (not equal)")
return 1
# convert a parallel vtu file into sequential one by glueing the pieces together
......@@ -365,7 +368,7 @@ def sort_vtk(root):
# sorts the data by point coordinates so that it is independent of index numbering
def sort_vtk_by_coordinates(root1, root2, verbose, convertedFromParallelVtu=False):
if not is_fuzzy_equal_node(root1.find(".//Points/DataArray"), root2.find(".//Points/DataArray"), absolute=1e-2, relative=1.5e-7, zeroValueThreshold=dict(), verbose=True, convertedFromParallelVtu=False):
if not is_fuzzy_equal_node(root1.find(".//Points/DataArray"), root2.find(".//Points/DataArray"), absolute=1e-2, relative=1.5e-7, zeroValueThreshold=dict(), verbose=False, convertedFromParallelVtu=False):
if verbose:
print("Sorting vtu by coordinates...")
for root in [root1, root2]:
......@@ -431,6 +434,40 @@ def sort_vtk_by_coordinates(root1, root2, verbose, convertedFromParallelVtu=Fals
for vertexIndex in cell:
largestCellMidPointForVertex[vertexIndex] = max(largestCellMidPointForVertex[vertexIndex], midpoint)
# floating point comparison operator for scalars
def float_cmp(a, b, eps):
if math.fabs(a-b) < eps:
return 0
elif a > b:
return 1
else:
return -1
# floating point comparison operator for vectors
def floatvec_cmp(a, b, eps):
for i, j in zip(a, b):
res = float_cmp(i, j, eps)
if res != 0:
return res
return 0
# compute an epsilon and a comparison operator for floating point comparisons
bBoxMax = max(vertexArray)
bBoxMin = min(vertexArray)
epsilon = math.sqrt(sum([(a-b)**2 for a, b in zip(bBoxMax, bBoxMin)]))*1e-7
# first compare by coordinates, if the same compare largestCellMidPointForVertex
# TODO: is there a more pythonic way?
def vertex_cmp(a, b):
res = floatvec_cmp(a[1], b[1], epsilon)
if res != 0:
return res
res2 = floatvec_cmp(largestCellMidPointForVertex[a[0]], largestCellMidPointForVertex[b[0]], epsilon)
if res2 != 0:
return res2
return 0
# obtain a vertex index map
vMap = []
for idx, coords in enumerate(vertexArray):
......@@ -439,7 +476,7 @@ def sort_vtk_by_coordinates(root1, root2, verbose, convertedFromParallelVtu=Fals
vertexIndexMap = [0]*len(vMap)
vertexIndexMapInverse = [0]*len(vMap)
# first sort by coordinates, if the same by largestCellMidPointForVertex
for idxNew, idxOld in enumerate(sorted(vMap, key=lambda x: (x[1], largestCellMidPointForVertex[x[0]]))):
for idxNew, idxOld in enumerate(sorted(vMap, key=functools.cmp_to_key(vertex_cmp))):
vertexIndexMap[idxOld[0]] = idxNew
vertexIndexMapInverse[idxNew] = idxOld[0]
......
......@@ -39,7 +39,7 @@ namespace Dumux {
/*!
* \ingroup InputOutput
* \brief Provides a grid manager with a method for creating creating vectors
* with polar Coordinates and one for creating a cartesian grid from
* with polar Coordinates and one for creating a Cartesian grid from
* these polar coordinates.
*/
template <class Grid>
......@@ -64,7 +64,7 @@ public:
const bool verbose = getParamFromGroup<bool>(modelParamGroup, "Grid.Verbosity", false);
std::array<std::vector<Scalar>, dim> polarCoordinates;
// Indices specifing in which direction the piece of cake is oriented
// Indices specifying in which direction the piece of cake is oriented
Dune::FieldVector<int, dim> indices(-1);
createVectors(polarCoordinates, indices, modelParamGroup, verbose);
......@@ -111,23 +111,24 @@ public:
if (hasRadial)
{
positions[i] = getParamFromGroup<std::vector<Scalar>>(modelParamGroup, "Grid.Radial" + std::to_string(i));
indices[0] = i; // Index specifing radial direction
indices[0] = i; // Index specifying radial direction
}
else if (hasAngular)
{
positions[i] = getParamFromGroup<std::vector<Scalar>>(modelParamGroup, "Grid.Angular" + std::to_string(i));
indices[1] = i; // Index specifing angular direction
indices[1] = i; // Index specifying angular direction
}
else // hasAxial
{
positions[i] = getParamFromGroup<std::vector<Scalar>>(modelParamGroup, "Grid.Axial" + std::to_string(i));
indices[2] = i; // Index specifing axial direction
indices[2] = i; // Index specifying axial direction
}
if (!std::is_sorted(positions[i].begin(), positions[i].end()))
DUNE_THROW(Dune::GridError, "Make sure to specify a monotone increasing \"Positions\" array");
if(positions[i].size() < 2)
DUNE_THROW(Dune::GridError, "Make sure to specify position arrays with at least two entries (min and max value).");
}
// check that all indices are specified i.e. > 0
......@@ -135,6 +136,14 @@ public:
if (indices[i] < 0)
DUNE_THROW(Dune::RangeError, "Please specify Positions Angular and Radial and Axial correctly and unambiguously!" << std::endl);
if(hasParam("Grid.WellRadius"))
{
std::cerr << "Deprecation warning: parameter Grid.WellRadius is deprecated. "
<< "Specify the WellRadius as the first radial coordinate." << std::endl;
positions[indices[0]][0] = getParamFromGroup<Scalar>(modelParamGroup, "Grid.WellRadius");
}
// the number of cells (has a default)
std::array<std::vector<int>, dim> cells;
for (int i = 0; i < dim; ++i)
......@@ -160,9 +169,10 @@ public:
using std::pow;
for (int dimIdx = 0; dimIdx < dim; dimIdx++)
{
std::size_t numGlobalPositions = 0;
// Each grid direction is subdivided into (numCells + 1) points
std::size_t numGlobalPositions = 1;
for (int zoneIdx = 0; zoneIdx < cells[dimIdx].size(); ++zoneIdx)
numGlobalPositions += cells[dimIdx][zoneIdx] + 1;
numGlobalPositions += cells[dimIdx][zoneIdx];
globalPositions[dimIdx].resize(numGlobalPositions);
std::size_t posIdx = 0;
......@@ -199,14 +209,15 @@ public:
increasingCellSize = true;
}
const bool useGrading = Dune::FloatCmp::eq(gradingFactor, 1.0) ? false : true;
// if the grading factor is exactly 1.0 do equal spacing
if (gradingFactor > 1.0 - 1e-7 && gradingFactor < 1.0 + 1e-7)
if (!useGrading)
{
height = 1.0 / numCells;
if (verbose)
std::cout << " -> h " << height * length << std::endl;
}
// if grading factor is not 1.0, do power law spacing
else
{
......@@ -226,7 +237,7 @@ public:
for (int i = 1; i < numCells; i++)
{
Scalar hI = height;
if (!(gradingFactor < 1.0 + 1e-7 && gradingFactor > 1.0 - 1e-7))
if (useGrading)
{
if (increasingCellSize)
hI *= pow(gradingFactor, i-1);
......@@ -259,10 +270,10 @@ public:
}
/*!
* \brief Creates cartesian grid from polar coordinates.
* \brief Creates Cartesian grid from polar coordinates.
*
* \param polarCoordinates Vector containing radial, angular and axial coordinates (in this order)
* \param indices Indices specifing the radial, angular and axial direction (in this order)
* \param indices Indices specifying the radial, angular and axial direction (in this order)
* \param modelParamGroup name of the model parameter group
* \param verbose if the output should be verbose
*/
......@@ -274,12 +285,12 @@ public:
std::vector<Scalar> dR = polarCoordinates[0];
std::vector<Scalar> dA = polarCoordinates[1];
// For the special case of 36o°, the last element is connected with the first one.
// For the special case of 360°, the last element is connected with the first one.
// Thus, one needs to distinguish here between all other cases, were the normal procedure
// is applied for all elements until the last one (= dA.size() - 1) and the 360° case,
// where it is only applied until the second to last element (dA.size() - 2).
int maxdA = dA.size() - 1;
if (Dune::FloatCmp::eq(polarCoordinates[1][polarCoordinates[1].size()-1], 360.0))
if (Dune::FloatCmp::eq(dA[dA.size()-1], 2*M_PI))
{
maxdA = dA.size() - 2;
}
......@@ -287,42 +298,60 @@ public:
GridFactory gridFactory;
constexpr auto type = Dune::GeometryTypes::cube(dim);
bool hasHole = true;
if(dR[0] < 1.0e-8*dR.back())
hasHole = false;
// create nodes
if (dim == 3)
{
constexpr auto prismType = Dune::GeometryTypes::prism;
std::vector<Scalar> dZ = polarCoordinates[2];
for (int j = 0; j <= dA.size() - 1; ++j)
for (int j = 0; j <= maxdA; ++j)
{
for (int l = 0; l <= dZ.size() - 1; ++l)
{
for (int i = 0; i <= dR.size()- 1; ++i)
for (int i = hasHole ? 0 : 1; i <= dR.size()- 1; ++i)
{
// Get radius for the well (= a hole) in the center
const auto wellRadius = getParamFromGroup<Scalar>(modelParamGroup, "Grid.WellRadius");
// transform into cartesian Coordinates
// transform into Cartesian coordinates
using std::cos;
using std::sin;
Dune::FieldVector <double, dim> v(0.0);
v[indices[2]] = dZ[l];
v[indices[0]] = cos(dA[j])*wellRadius + cos(dA[j])*dR[i];
v[indices[1]] = sin(dA[j])*wellRadius + sin(dA[j])*dR[i];
v[indices[0]] = cos(dA[j])*dR[i];
v[indices[1]] = sin(dA[j])*dR[i];
if (verbose)
{
std::cout << "Coordinates of : "
<< v[0] << " " << v[1] << " " << v[2] << std::endl;
}
printCoordinate(v);
gridFactory.insertVertex(v);
}
}
}
// if the well is not considered, we have to add the points lying on the center-line
// this has to be done separately, otherwise these points would occur multiple times
if(!hasHole)
{
for (int l = 0; l <= dZ.size() - 1; ++l)
{
Dune::FieldVector <double, dim> v(0.0);
v[indices[2]] = dZ[l];
v[indices[0]] = 0.0;
v[indices[1]] = 0.0;
if (verbose)
printCoordinate(v);
gridFactory.insertVertex(v);
}
}
if (verbose)
std::cout << "Filled node vector" << std::endl;
// assign nodes
unsigned int z = 0;
unsigned int t = 0;
unsigned int rSize = hasHole ? dR.size() : dR.size()-1;
unsigned int zSize = dZ.size();
for (int j = 0; j < dA.size() - 1; ++j)
{
for (int l = 0; l < dZ.size() - 1; ++l)
......@@ -331,130 +360,171 @@ public:
{
for (int i = 0; i < dR.size() - 1; ++i)
{
unsigned int rSize = dR.size();
unsigned int zSize = dZ.size();
std::vector<unsigned int> vid({z, z+1, z+rSize*zSize,
z+rSize*zSize+1, z+rSize, z+rSize+1,
z+rSize*zSize+rSize, z+rSize*zSize+rSize+1});
if (verbose)
if(!hasHole && i==0)
{
std::cout << "element vertex ids: ";
for (int k = 0; k < vid.size(); ++k)
std::cout << vid[k] << " ";
std::cout << std::endl;
std::vector<unsigned int> vid({rSize*zSize*(maxdA+1) + l, z, z+rSize*zSize,
rSize*zSize*(maxdA+1) + l+1, z+rSize, z+rSize*zSize+rSize});
if (verbose)
printIndices(vid);
gridFactory.insertElement(prismType, vid);
}
else
{
std::vector<unsigned int> vid({z, z+1,
z+rSize*zSize, z+rSize*zSize+1,
z+rSize, z+rSize+1,
z+rSize*zSize+rSize, z+rSize*zSize+rSize+1});
gridFactory.insertElement(type, vid);
if (verbose)
printIndices(vid);
gridFactory.insertElement(type, vid);
z = z+1;
z++;
}
}
z = z+1;
z++;
}
else
{
// assign nodes for 360°-cake
for (int i = 0; i < dR.size() - 1; ++i)
{
// z = z + 1;
unsigned int rSize = dR.size();
std::vector<unsigned int> vid({z, z+1, t,
t+1, z+rSize, z+rSize+1,
t+rSize, t+rSize+1});
if (verbose)
if(!hasHole && i==0)
{
std::cout << "element vertex ids: ";
for (int k = 0; k < vid.size(); ++k)
std::cout << vid[k] << " ";
std::cout << std::endl;
std::vector<unsigned int> vid({rSize*zSize*(maxdA+1) + l, z, t,
rSize*zSize*(maxdA+1) + l+1, z+rSize, t+rSize});
if (verbose)
printIndices(vid);
gridFactory.insertElement(prismType, vid);
}
else
{
std::vector<unsigned int> vid({z, z+1,
t, t+1,
z+rSize, z+rSize+1,
t+rSize, t+rSize+1});
gridFactory.insertElement(type, vid);
t = t + 1;
z = z+1;
if (verbose)
printIndices(vid);
gridFactory.insertElement(type, vid);
t++;
z++;
}
}
t = t + 1;
z = z+1;
t++;
z++;
if (verbose)
std::cout << "assign nodes 360° ends..." << std::endl;
}
if (verbose)
std::cout << "assign nodes 360° ends..." << std::endl;
}
z = z + dR.size();
z += rSize;
}
}
// for dim = 2
else
{
for (int j = 0; j <= dA.size() - 1; ++j)
constexpr auto triangleType = Dune::GeometryTypes::simplex(dim);
for (int j = 0; j <= maxdA; ++j)
{
for (int i = 0; i <= dR.size()- 1; ++i)
for (int i = hasHole ? 0 : 1; i <= dR.size()- 1; ++i)
{
// Get radius for the well (= a hole) in the center
const Scalar wellRadius = getParamFromGroup<Scalar>(modelParamGroup, "Grid.WellRadius");
// transform into cartesian Coordinates
// transform into Cartesian coordinates
Dune::FieldVector <double, dim> v(0.0);
v[indices[0]] = cos(dA[j])*wellRadius + cos(dA[j])*dR[i];
v[indices[1]] = sin(dA[j])*wellRadius + sin(dA[j])*dR[i];
if(verbose) std::cout << "Coordinates of : " << v[0] << " " << v[1] << std::endl;
v[indices[0]] = cos(dA[j])*dR[i];
v[indices[1]] = sin(dA[j])*dR[i];
if(verbose)
printCoordinate(v);
gridFactory.insertVertex(v);
}
}
// if the well is not considered, we have to add the point located at the center
// this has to be done separately, otherwise this point would occur multiple times
if(!hasHole)
{
Dune::FieldVector <double, dim> v(0.0);
v[indices[0]] = 0.0;
v[indices[1]] = 0.0;
if(verbose)
printCoordinate(v);
gridFactory.insertVertex(v);
}
std::cout << "Filled node vector" << std::endl;
// assign nodes
unsigned int z = 0;
unsigned int t = 0;
unsigned int rSize = hasHole ? dR.size() : dR.size()-1;
for (int j = 0; j < dA.size() - 1; ++j)
{
if (j < maxdA)
{
for (int i = 0; i < dR.size() - 1; ++i)
{
unsigned int rSize = dR.size();
std::vector<unsigned int> vid({z, z+1, z+rSize, z+rSize+1});
if (verbose)
if(!hasHole && i==0)
{
std::cout << "element vertex ids: ";
for (int k = 0; k < vid.size(); ++k)
std::cout << vid[k] << " ";
std::cout << std::endl;
std::vector<unsigned int> vid({rSize*(maxdA+1), z, z+rSize});
if (verbose)
printIndices(vid);
gridFactory.insertElement(triangleType, vid);
}
else
{
std::vector<unsigned int> vid({z, z+1, z+rSize, z+rSize+1});
gridFactory.insertElement(type, vid);
z = z+1;
if (verbose)
printIndices(vid);
gridFactory.insertElement(type, vid);
z++;
}
}
z = z+1;
z++;
}
else
{
// assign nodes for 360°-cake
for (int i = 0; i < dR.size() - 1; ++i)
{
// z = z + 1;
std::vector<unsigned int> vid({z, z+1, t, t+1});
if (verbose)
if(!hasHole && i==0)
{
std::cout << "element vertex ids: ";
for (int k = 0; k < vid.size(); ++k)
std::cout << vid[k] << " ";
std::cout << std::endl;
std::vector<unsigned int> vid({rSize*(maxdA+1), z, t});
if (verbose)
printIndices(vid);
gridFactory.insertElement(triangleType, vid);
}
else
{
std::vector<unsigned int> vid({z, z+1, t, t+1});
if (verbose)
printIndices(vid);
gridFactory.insertElement(type, vid);
t = t + 1;
z = z+1;
gridFactory.insertElement(type, vid);
t++;
z++;
}
}
t = t + 1;
z = z+1;
t++;
z++;
if (verbose)
std::cout << "assign nodes 360 ends..." << std::endl;
}
std::cout << "assign nodes 360 ends..." << std::endl;
}
}
// return the grid pointer
......@@ -479,6 +549,21 @@ public:
}
protected:
static void printCoordinate(const Dune::FieldVector <double, dim>& v)
{
std::cout << "Coordinates of : ";
for (int k = 0; k < v.size(); ++k)
std::cout << v[k] << " ";
std::cout << std::endl;
}
static void printIndices(const std::vector<unsigned int>& vid)
{
std::cout << "element vertex indices: ";
for (int k = 0; k < vid.size(); ++k)
std::cout << vid[k] << " ";
std::cout << std::endl;
}
/*!
* \brief Returns a reference to the shared pointer to the grid.
......
add_input_file_links()
dune_symlink_to_source_files(FILES grids)
add_executable(test_gridmanager_cake EXCLUDE_FROM_ALL test_gridmanager_cake.cc)
add_executable(test_gridmanager_cake_ug EXCLUDE_FROM_ALL test_gridmanager_cake.cc)
target_compile_definitions(test_gridmanager_cake_ug PUBLIC "USEUG=1")
dune_add_test(NAME test_gridmanager_cake_360
TARGET test_gridmanager_cake
add_executable(test_gridmanager_cake_alu EXCLUDE_FROM_ALL test_gridmanager_cake.cc)
dune_add_test(NAME test_gridmanager_cake_360_ug
TARGET test_gridmanager_cake_ug
LABELS unit io
CMAKE_GUARD "( dune-uggrid_FOUND OR dune-alugrid_FOUND )"
CMAKE_GUARD dune-uggrid_FOUND
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--command "${CMAKE_CURRENT_BINARY_DIR}/test_gridmanager_cake -Grid.Name 360"
--command "${CMAKE_CURRENT_BINARY_DIR}/test_gridmanager_cake_ug -Grid.Name ug-360"
--files ${CMAKE_SOURCE_DIR}/test/references/test_gridmanager_cake_3d_360-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/cake-3d-360.vtu
${CMAKE_CURRENT_BINARY_DIR}/cake-3d-ug-360.vtu
${CMAKE_SOURCE_DIR}/test/references/test_gridmanager_cake_2d_360-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/cake-2d-360.vtu)
${CMAKE_CURRENT_BINARY_DIR}/cake-2d-ug-360.vtu
--relative 1e-6)
dune_add_test(NAME test_gridmanager_cake_360_alu
TARGET test_gridmanager_cake_alu
LABELS unit io
CMAKE_GUARD dune-alugrid_FOUND
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--command "${CMAKE_CURRENT_BINARY_DIR}/test_gridmanager_cake_alu -Grid.Name alu-360"
--files ${CMAKE_SOURCE_DIR}/test/references/test_gridmanager_cake_3d_360-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/cake-3d-alu-360.vtu
${CMAKE_SOURCE_DIR}/test/references/test_gridmanager_cake_2d_360-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/cake-2d-alu-360.vtu
--relative 1e-6)
dune_add_test(NAME test_gridmanager_cake_210_ug
TARGET test_gridmanager_cake_ug
LABELS unit io
CMAKE_GUARD dune-uggrid_FOUND
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--command "${CMAKE_CURRENT_BINARY_DIR}/test_gridmanager_cake_ug -Grid.Name ug-210 -Grid.Angular1 '0.0 210.0'"
--files ${CMAKE_SOURCE_DIR}/test/references/test_gridmanager_cake_3d_210-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/cake-3d-ug-210.vtu
${CMAKE_SOURCE_DIR}/test/references/test_gridmanager_cake_2d_210-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/cake-2d-ug