Commit 96fa0112 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

cleanup of gridfunctions

parent 75106fc1
...@@ -4,7 +4,6 @@ import collections ...@@ -4,7 +4,6 @@ import collections
x, y, z = symbols("x y z") x, y, z = symbols("x y z")
a, b, c = symbols("a b c") a, b, c = symbols("a b c")
u, v = symbols("u v")
# simplify entries of an array # simplify entries of an array
def simplify_all(A): def simplify_all(A):
...@@ -50,19 +49,17 @@ v = acos(X0[2]) ...@@ -50,19 +49,17 @@ v = acos(X0[2])
X1 = [a*cos(u)*sin(v), b*sin(u)*sin(v), c*cos(v)] X1 = [a*cos(u)*sin(v), b*sin(u)*sin(v), c*cos(v)]
# jacobian of parametrization # jacobian of parametrization
J = simplify_all([[diff(X1[i],X[j]) for i in range(3)] for j in range(3)]) J = [[diff(X1[i],X[j]) for i in range(3)] for j in range(3)]
print("J = ") print("J = ")
print("return {") print("return {")
for i in range(3): for i in range(3):
print(" {") print(" {")
for j in range(3): for j in range(3):
print(" ",ccode(J[i][j]),",") print(" ",ccode(simplify(J[i][j])),",")
print(" },") print(" },")
print("};") print("};")
print("") print("")
print("equal = ",simplify((x**2 + y**2 + z**2)**2 - (x**4 + 2*x**2*y**2 + 2*x**2*z**2 + y**4 + 2*y**2*z**2 + z**4)))
# P(F)_j # P(F)_j
def project1(F): def project1(F):
return simplify_all([ return simplify_all([
......
...@@ -5,9 +5,11 @@ ...@@ -5,9 +5,11 @@
#include <type_traits> #include <type_traits>
#include <dune/curvedsurfacegrid/gridfunctions/analyticgridfunction.hh> #include <dune/common/math.hh>
#include <dune/functions/common/defaultderivativetraits.hh> #include <dune/functions/common/defaultderivativetraits.hh>
#include "analyticgridfunction.hh"
namespace Dune namespace Dune
{ {
// Ellipsoid functor // Ellipsoid functor
...@@ -27,6 +29,7 @@ namespace Dune ...@@ -27,6 +29,7 @@ namespace Dune
{} {}
//! project the coordinate to the ellipsoid //! project the coordinate to the ellipsoid
// NOTE: This is not a closes-point projection, but a spherical-coordinate projection
template< class Domain > template< class Domain >
Domain operator() (const Domain& X) const Domain operator() (const Domain& X) const
{ {
...@@ -36,11 +39,11 @@ namespace Dune ...@@ -36,11 +39,11 @@ namespace Dune
} }
//! derivative of the projection //! derivative of the projection
#if0 friend auto derivative (const EllipsoidProjection& ellipsoid)
friend auto derivative (const EllipsoidProjection& sphere)
{ {
return [r=sphere.radius_](auto const& X) return [a=ellipsoid.a_,b=ellipsoid.b_,c=ellipsoid.c_](auto const& X)
{ {
using std::sqrt;
using Domain = std::decay_t<decltype(X)>; using Domain = std::decay_t<decltype(X)>;
using DerivativeTraits = Functions::DefaultDerivativeTraits<Domain(Domain)>; using DerivativeTraits = Functions::DefaultDerivativeTraits<Domain(Domain)>;
typename DerivativeTraits::Range out; typename DerivativeTraits::Range out;
...@@ -53,8 +56,8 @@ namespace Dune ...@@ -53,8 +56,8 @@ namespace Dune
T nrm1 = x2 + y2 + z2; T nrm1 = x2 + y2 + z2;
T nrm2 = sqrt(nrm0/nrm1); T nrm2 = sqrt(nrm0/nrm1);
T nrm3 = sqrt(nrm0/x2); T nrm3 = sqrt(nrm0/x2);
T nrm4 = power(sqrt(nrm1), 3); T nrm4 = sqrt(nrm1)*nrm1;
T nrm5 = sqrt(nrm0)*nrm1*nrm1/sqrt(nrm1); T nrm5 = sqrt(nrm0)*nrm4;
return { return {
{ {
...@@ -73,11 +76,8 @@ namespace Dune ...@@ -73,11 +76,8 @@ namespace Dune
c*(x2 + y2)/nrm4 c*(x2 + y2)/nrm4
} }
}; };
}; };
} }
#endif
//! Normal vector //! Normal vector
template< class Domain > template< class Domain >
......
...@@ -7,17 +7,17 @@ namespace Dune ...@@ -7,17 +7,17 @@ namespace Dune
{ {
//! A set of entities of given `codim` of a `Grid` //! A set of entities of given `codim` of a `Grid`
/** /**
* \tparam G The grid type * \tparam GridType The grid type
* \tparam codim Codimension of the entities to define the set of. * \tparam codim Codimension of the entities to define the set of.
* *
* \note This entityset just defines types * \note This entityset just defines types
**/ **/
template< class G, int codim > template< class GridType, int codim >
class GridEntitySet class GridEntitySet
{ {
public: public:
//! Type of the grid //! Type of the grid
using Grid = G; using Grid = GridType;
//! Type of Elements contained in this EntitySet //! Type of Elements contained in this EntitySet
using Element = typename Grid::template Codim<codim>::Entity; using Element = typename Grid::template Codim<codim>::Entity;
......
...@@ -28,14 +28,18 @@ namespace Dune ...@@ -28,14 +28,18 @@ namespace Dune
} // end namespace Impl } // end namespace Impl
//! Type-Traits to extract the grid-type from a GridFunction `GF`
template< class GF > template< class GF >
using GridOf_t = typename Impl::GridOf<std::decay_t<decltype(std::declval<GF>().entitySet())>>::type; using GridOf_t = typename Impl::GridOf<std::decay_t<decltype(std::declval<GF>().entitySet())>>::type;
//! Conditionally define the GridFunction type.
template< class HostGrid, class GF > /**
* If the type `GF` is a GridFunction on the grid `Grid`, use this type as GridFunction type,
* Otherwise it is a functor and construct an \ref AnalyticGridFunction on the `Grid`.
**/
template< class Grid, class GF >
using GridFunctionOf_t using GridFunctionOf_t
= std::conditional_t<Dune::Concept::isGridFunction<GF, HostGrid>(), GF, AnalyticGridFunction<HostGrid,GF>>; = std::conditional_t<Dune::Concept::isGridFunction<GF, Grid>(), GF, AnalyticGridFunction<Grid,GF>>;
} // end namespace Dune } // end namespace Dune
......
...@@ -75,7 +75,7 @@ namespace Dune ...@@ -75,7 +75,7 @@ namespace Dune
} }
T mean_curvature (FieldVector<T,3> X) T mean_curvature (FieldVector<T,3> X) const
{ {
using std::sqrt; using std::sqrt;
X = (*this)(X); X = (*this)(X);
......
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