#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 dune-grid 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 i-th 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 cross-product 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==cdim-1 -
GlobalCoordinate integrateAnalyticalTimes(const PolynomialVector & PVec)
- integrates (PVec x normal) over the surface, when mydim==cdim-1
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)