diff --git a/dune/fufem/test/dunepythontest.cc b/dune/fufem/test/dunepythontest.cc index 34abf538f0c0b8a926780e17772a6c9fef16bf3f..86798acf8b7914d32c9c8b782a754f388e878824 100644 --- a/dune/fufem/test/dunepythontest.cc +++ b/dune/fufem/test/dunepythontest.cc @@ -13,6 +13,13 @@ #include <dune/common/function.hh> #include <dune/common/parametertree.hh> +#if HAVE_UG +#include <dune/grid/uggrid.hh> +#include <dune/grid/uggrid/uggridfactory.hh> +#include <dune/grid/common/gridfactory.hh> +#include <dune/grid/io/file/vtk/vtkwriter.hh> +#endif + #include <dune/functions/common/differentiablefunctionfromcallables.hh> #include <dune/fufem/dunepython.hh> @@ -759,6 +766,35 @@ int main(int argc, char** argv) std::cout << "Result of calling a function double(double) : " << h(0.5) << std::endl; } + +#if HAVE_UG + std::cout << std::endl; + std::cout << "Example 17: Filling a GridFactory form python **********************************" << std::endl; + { + // prepare types and grid factory + using Grid = Dune::UGGrid<2>; + using GridFactory = Dune::GridFactory<Grid>; + GridFactory gridFactory; + + // Get grid describtion from python and convert it to GridFactory + auto pyGrid = pyMain.get("grid"); + pyGrid.toC(gridFactory); + + // Let GridFactory create grid + std::unique_ptr<Grid> gridPtr(gridFactory.createGrid()); + auto& grid = *gridPtr; + + // Refine and write grid to visualize parametrized boundary + grid.globalRefine(1); + grid.globalRefine(1); + + Dune::VTKWriter<Grid::LeafGridView> vtkWriter(grid.leafGridView()); + vtkWriter.write("python-generated-grid"); + + std::cout << "number of elements in created grid : " << grid.leafGridView().indexSet().size(0) << std::endl; + } +#endif + // Before exiting you should stop the python interpreter. Python::stop(); diff --git a/dune/fufem/test/dunepythontest.py b/dune/fufem/test/dunepythontest.py index bebcc3cc8fa86c2b3dcdeb4bbef7df1f6c7a9db2..c82d9397b04326b03732daeca5d13c8c52adfebb 100644 --- a/dune/fufem/test/dunepythontest.py +++ b/dune/fufem/test/dunepythontest.py @@ -1,3 +1,6 @@ +from numpy import * +from numpy.linalg import norm + foo='foo' bar='bar' @@ -9,3 +12,78 @@ str_list=["string1","string2","string3"] int_tuple=(1,2,3) + +# some utilities for creating parametrized boundaries + +class BoundarySegment(list): + """ A python BoundarySegment representation + + It behaves like a list and function at the same time. + You can access the passed vertices using [] and + call the passed parametrization usin (). + """ + def __init__(self, vertices, parametrization): + super(BoundarySegment,self).__init__(vertices) + self.parametrization = parametrization + + def __call__(self, x): + return self.parametrization(x) + +class TransformedBoundarySegment(BoundarySegment): + """ This transforms a local parametrization into global one + + The local parametrization is assumed to be a curve f + on [0,1] with f(0)=(0,0) and f(1)=(1,0). + This class rotatetes, scales, and translates the graph + of the curve such that the transformed curve Tf satisfies + Tf(0) = coordinates(vertices[0]) and + Tf(1) = coordinates(vertices[1]). + The scaling in normal direction is done with the + same factor such that T is angle and curvature preserving. + """ + def __init__(self, vertices, localParametrization, allCoordinates): + coordinates = array([allCoordinates[i] for i in vertices]) + edge = coordinates[1]-coordinates[0] + normal = (-edge[1], edge[0]) + parametrization = lambda x: dot(transpose(array((edge, normal))), localParametrization(x)) + coordinates[0] + super(TransformedBoundarySegment, self).__init__(vertices, parametrization) + +class ArcOverSecant: + """ Parametrization of a circle arc over a secant + + The circle arc is parametrized over [0,1] and + placed such that its secant line goes from + (0,0) to (1,0). The radius of the arc is + selected automatically such that the + arc has the given angle. + """ + def __init__(self, angle): + self.alphaHalf = float(angle)/2.0 + self.invSinAlphaHalf = 1.0/sin(self.alphaHalf) + + def __call__(self, x): + betaHalf = self.alphaHalf*x[0] + secantlength = sin(betaHalf)*self.invSinAlphaHalf + secantangle = self.alphaHalf-betaHalf + return (cos(secantangle)*secantlength, sin(secantangle)*secantlength) + + + + + +# reentrent corner circle grid + +class Grid: + def __init__(self): + self.vertices = [(0,0), (1,0), (0,1), (-1,0), (0, -1)] + self.elements = [(0,1,2), (0,2,3), (0,3,4)] + + # We use some utility classes here, but each segment essentially needs to be iterable, + # e.g, a list. If it is additionally callable, it is used as parametrization. + self.segments = [] + self.segments.append(TransformedBoundarySegment((2,1), ArcOverSecant(2.0*pi/4), self.vertices)) + self.segments.append(TransformedBoundarySegment((3,2), ArcOverSecant(2.0*pi/4), self.vertices)) + self.segments.append(TransformedBoundarySegment((4,3), ArcOverSecant(2.0*pi/4), self.vertices)) + +grid = Grid() +