#1583 Additional user grid.geometry() methods from CurvilinearGeometry
Metadata
Property  Value 

Reported by  Aleksejs Fomins (aleksejs.fomins@lspr.ch) 
Reported at  Mar 5, 2015 12:00 
Type  Feature Request 
Version  2.3 
Operating System  Unspecified / All 
Description
In this post I will discuss some of the new functionality available in CurvilinearGeometry
. The question is to find the best way to provide the user with this functionality through dunegrid interface. This post will be quite big, I foresee to split it into smaller requests at a later stage
Priority 1
methods to obtain basic information about curvilinear entity

InterpolationOrder()
 Gets polynomial order of this entity 
vertex(int i)
 Gets ith interpolatory vertex coordinate (likecorner(int i)
which provides corners of the geometry) 
vertices()
 Gets number of interpolatory vertices (likecorners()
which provides the number of corners) 
isInside(local, global)
 already mentioned in request #1578
Priority 2
When investigating the available quadrature rules, we did not find a suitable way to integrate polynomials of high order over triangles and tetrahedra
To go around this problem, we implemented a class which allows to analytically represent polynomials, and analytically integrate them over the reference elements. That way an integral of a polynomial function over a curvilinear entity can be done analytically for arbitrary order as long as the system has enough resources. This approach has 2 downsides and 1 additional upside
 (Currently) it only works if the integrand is polynomial. Further investigation is necessary if some other basis is used as integrand
 It does not work when
mydim < cdim
. When integrating a scalar function over, e.g., a surface in 3D, the Jacobian determinant is a square root of a polynomial, so in this case numerical approach is necessary.  It is, however, possible, to find surface integrals (dot and crossproduct with normal) analytically.
so the available methods are

ctype integrateAnalyticalScalar(const LocalPolynomial & P)
 integrates a polynpmial over the entity when mydim==cdim 
ctype integrateAnalyticalDot(const PolynomialVector & PVec)
 integrates (PVec.normal) over the surface, when mydim==cdim1 
GlobalCoordinate integrateAnalyticalTimes(const PolynomialVector & PVec)
 integrates (PVec x normal) over the surface, when mydim==cdim1
Since we are interested to be able to integrate polynomial scalars over all codimensions 3D, we provide an internal iterative scheme to do so over triangles and edges. Our scheme can integrate arbitrary functor f()(x) over the curvilinear triangle or edge (aka integrating f(x)sqrt(det(J(x))) over reference edge/triangle). It works by recursively refining the entity, integrating a 2nd and 4th order interpolating order polynomial over each leaf, and taking the difference of the two as the running error estimate. We have tested it and it is quite precise, but rather slow for high orders. Any contributions welcome
So the numerical methods are

volume()
 standard dune.geometry method requires numerical integration in case mydim < cdim, which is done internally 
integrateNumerical(Functor f)
 iterative scheme to integrate arbitrary functor over simplex of dimension 1 or 2 
integrateNumerical
needs to be extended to manage dim=3 arbitrary integrals, it is not available at the moment, since I was not sure if it is worth putting extra effort into the method since it is not the best algorithm there is anyway
Priority 3  auxiliary polynomial functions
There are some auxiliary functions that may or may not be interesting to the user

PolynomialVector interpolatoryVectorAnalytical()
 Analytical representation of the curvilinear local>global map using a vector of polynomials of dimension cdim 
LocalPolynomial JacobianDeterminantAnalytical()
 Analytical representation of Jacobian determinant when mydim=cdim 
PolynomialVector NormalIntegrationElementAnalytical()
 vec{dS} in 2D 
LocalPolynomial IntegrationElementSquaredAnalytical()
 det(J^T J)