Commit 2133b97f authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

grid implementation cleaned up

parent 18d9ad99
......@@ -47,6 +47,12 @@ X0 = [x/nrm_X, y/nrm_X, z/nrm_X]
u = atan(X0[1]/X0[0])
v = acos(X0[2])
X1 = [a*cos(u)*sin(v), b*sin(u)*sin(v), c*cos(v)]
print("X1 = ")
print("return {")
for i in range(3):
print(" ",ccode(X1[i]),",")
print("};")
print("")
# jacobian of parametrization
J = [[diff(X1[i],X[j]) for i in range(3)] for j in range(3)]
......
......@@ -29,68 +29,64 @@ namespace Dune
struct hasSingleGeometryType< CurvedSurfaceGrid<GridFunction,order> >
{
using HostGrid = GridOf_t<GridFunction>;
static const bool v = hasSingleGeometryType< HostGrid > :: v;
static const unsigned int topologyId = hasSingleGeometryType< HostGrid > :: topologyId;
static const bool v = hasSingleGeometryType<HostGrid>::v;
static const unsigned int topologyId = hasSingleGeometryType<HostGrid>::topologyId;
};
template< class GridFunction, int order, int codim >
struct hasEntity< CurvedSurfaceGrid<GridFunction,order>, codim >
{
static const bool v = true;
using HostGrid = GridOf_t<GridFunction>;
static const bool v = hasEntity<HostGrid,codim>::v;
};
template< class GridFunction, int order, int codim >
struct hasEntityIterator< CurvedSurfaceGrid<GridFunction,order>, codim >
{
static const bool v = true;
using HostGrid = GridOf_t<GridFunction>;
static const bool v = hasEntityIterator<HostGrid,codim>::v;
};
//! Implements geometry only for codim=0 entity
template< class GridFunction, int order, int codim >
struct canCommunicate< CurvedSurfaceGrid<GridFunction,order>, codim >
struct hasGeometry< CurvedSurfaceGrid<GridFunction,order>, codim >
{
using HostGrid = GridOf_t<GridFunction>;
static const bool v = canCommunicate< HostGrid, codim >::v && hasEntity< HostGrid, codim >::v;
static const bool v = (codim == 0);
};
template< class GridFunction, int order >
struct hasBackupRestoreFacilities< CurvedSurfaceGrid<GridFunction,order> >
template< class GridFunction, int order, int codim >
struct canCommunicate< CurvedSurfaceGrid<GridFunction,order>, codim >
{
using HostGrid = GridOf_t<GridFunction>;
static const bool v = hasBackupRestoreFacilities< HostGrid >::v;
static const bool v = canCommunicate<HostGrid, codim>::v && hasEntity<HostGrid, codim>::v;
};
//! Conformity of the grid is not guaranteed since it depends on the GridFunction that
//! must be continuous in that case.
template< class GridFunction, int order >
struct isLevelwiseConforming< CurvedSurfaceGrid<GridFunction,order> >
{
using HostGrid = GridOf_t<GridFunction>;
static const bool v = isLevelwiseConforming< HostGrid >::v;
static const bool v = false;
};
//! Conformity of the grid is not guaranteed since it depends on the GridFunction that
//! must be continuous in that case.
template< class GridFunction, int order >
struct isLeafwiseConforming< CurvedSurfaceGrid<GridFunction,order> >
{
using HostGrid = GridOf_t<GridFunction>;
static const bool v = isLeafwiseConforming< HostGrid >::v;
};
template< class GridFunction, int order >
struct threadSafe< CurvedSurfaceGrid<GridFunction,order> >
{
static const bool v = false;
};
//! Implements only partial backup-restore facilities, since the gridfunction is not
//! backuped automatically.
template< class GridFunction, int order >
struct viewThreadSafe< CurvedSurfaceGrid<GridFunction,order> >
struct hasBackupRestoreFacilities< CurvedSurfaceGrid<GridFunction,order> >
{
static const bool v = false;
using HostGrid = GridOf_t<GridFunction>;
static const bool v = hasBackupRestoreFacilities<HostGrid>::v;
};
// hasHostEntity
// -------------
......@@ -98,7 +94,7 @@ namespace Dune
struct hasHostEntity< CurvedSurfaceGrid<GridFunction,order>, codim >
{
using HostGrid = GridOf_t<GridFunction>;
static const bool v = hasEntity< HostGrid, codim >::v;
static const bool v = hasEntity<HostGrid, codim>::v;
};
} // namespace Capabilities
......
......@@ -4,53 +4,70 @@
#define DUNE_CURVED_SURFACE_GRID_CONCEPTS_HH
#include <dune/common/concept.hh>
#include <dune/curvedsurfacegrid/gridfunctions/gridentityset.hh>
#include <dune/functions/common/functionconcepts.hh>
namespace Dune {
namespace Concept {
//! Concept: objects that can be called with given argument list `Args...`
template< class... Args >
struct Callable
{
template< class F >
auto require (F&& f) -> decltype(
f(std::declval<Args>()...)
);
};
//! Check if `F` models the `Callable` concept with the given arguments `Args...`
template< class F, class... Args >
constexpr auto isCallable ()
{ return models<Concept::Callable<Args...>, F>(); }
//! Concept: function associated to the `LocalContext` that can be evaluated in local coordinates
template< class LocalContext >
struct LocalFunction
{
using LocalCoordinate = typename LocalContext::Geometry::LocalCoordinate;
template< class LF >
auto require(LF&& lf) -> decltype(
auto require (LF&& lf) -> decltype(
lf.bind(std::declval<LocalContext>()),
lf.unbind(),
lf.localContext(),
requireConcept<Dune::Functions::Concept::Callable<LocalCoordinate>>(lf),
requireConcept<Concept::Callable<LocalCoordinate>>(lf),
requireConvertible<LocalContext>(lf.localContext())
);
};
//! Check if `LF` models the `LocalFunction` concept on the given `LocalContext`
template< class LF, class LocalContext >
constexpr bool isLocalFunction()
constexpr bool isLocalFunction ()
{ return models<Concept::LocalFunction<LocalContext>, LF>(); }
template< class HostGrid >
//! Concept: function associated to the `Grid` that can be evaluated in global coordinates
template< class Grid >
struct GridFunction
{
using EntitySet = GridEntitySet<HostGrid,0>;
using LocalContext = typename EntitySet::Element;
using GlobalCoordinate = typename EntitySet::GlobalCoordinate;
using LocalContext = typename Grid::template Codim<0>::Entity;
using GlobalCoordinate = typename LocalContext::Geometry::GlobalCoordinate;
template< class GF >
auto require(GF&& gf) -> decltype(
auto require (GF&& gf) -> decltype(
localFunction(gf),
gf.entitySet(),
requireConcept<Dune::Functions::Concept::Callable<GlobalCoordinate>>(gf),
requireConcept<LocalFunction<LocalContext>>(localFunction(gf))
requireConcept<Concept::Callable<GlobalCoordinate>>(gf),
requireConcept<Concept::LocalFunction<LocalContext>>(localFunction(gf))
);
};
template< class GF, class HostGrid >
constexpr bool isGridFunction()
{ return models<Concept::GridFunction<HostGrid>, GF>(); }
//! Check if `GF` models the `GridFunction` concept on the given `Grid`
template< class GF, class Grid >
constexpr bool isGridFunction ()
{ return models<Concept::GridFunction<Grid>, GF>(); }
} // end namespace Concept
} // end namespace Dune
} // end namespace Dune
#endif // DUNE_CURVED_SURFACE_GRID_CONCEPTS_HH
\ No newline at end of file
......@@ -5,6 +5,4 @@
#include <dune/curvedsurfacegrid/grid.hh>
#include <dune/grid/geometrygrid/persistentcontainer.hh>
// add your classes here
#endif // DUNE_CURVEDSURFACEGRID_HH
......@@ -14,7 +14,6 @@
#include <dune/curvedsurfacegrid/datahandle.hh>
#include <dune/curvedsurfacegrid/gridfunctions/analyticgridfunction.hh>
#include <dune/curvedsurfacegrid/gridfunctions/gridfunction.hh>
#include <dune/functions/common/functionconcepts.hh>
#include <dune/grid/geometrygrid/identity.hh>
#include <dune/grid/geometrygrid/persistentcontainer.hh>
#include <dune/grid/geometrygrid/grid.hh>
......@@ -59,21 +58,21 @@ namespace Dune
class CurvedSurfaceGrid;
//! Generator for CurvedSurfaceGrid from a grid-functions
template< class HostGrid, class GF, int ORDER = -1,
std::enable_if_t<Dune::Concept::isGridFunction<GF, HostGrid>(), int> = 0>
auto curedSurfaceGrid (HostGrid& hostGrid, GF&& gridFunction)
template< class HostGrid, class GF, int ORDER = -1,
std::enable_if_t<Concept::isGridFunction<GF, HostGrid>(), int> = 0>
auto curedSurfaceGrid (HostGrid& hostGrid, GF&& gridFunction)
{
static_assert(std::is_same<HostGrid, GridOf_t<std::decay_t<GF>>>::value, "GridFunction must be defined on the HostGrid");
return CurvedSurfaceGrid<std::decay_t<GF>,ORDER>{hostGrid, std::forward<GF>(gridFunction)};
}
//! Generator for CurvedSurfaceGrid from a callable
template< class HostGrid, class F, int ORDER = -1,
std::enable_if_t<not Dune::Concept::isGridFunction<F, HostGrid>(), int> = 0>
auto curedSurfaceGrid (HostGrid& hostGrid, F&& callable)
template< class HostGrid, class F, int ORDER = -1,
std::enable_if_t<not Concept::isGridFunction<F, HostGrid>(), int> = 0>
auto curedSurfaceGrid (HostGrid& hostGrid, F&& callable)
{
using GlobalCoordinate = typename GridEntitySet<HostGrid,0>::GlobalCoordinate;
static_assert(Functions::Concept::isCallable<F, GlobalCoordinate>(), "Function must be callable");
static_assert(Concept::isCallable<F, GlobalCoordinate>(), "Function must be callable");
auto gridFct = analyticGridFunction<HostGrid>(std::forward<F>(callable));
return CurvedSurfaceGrid<decltype(gridFct),ORDER>{hostGrid, std::move(gridFct)};
}
......@@ -217,7 +216,7 @@ namespace Dune
* \param[in] hostGrid reference to the grid to wrap
* \param[in] gridFunction mapping from global coordinates in the host geometry
* to global coordinates in the curved geometry
*
*
* If the gridFunction is passed by (non-const) reference it is stored as a non-destroying shared_ptr.
* Otherwise it is copied or moved into a new object (stored as shared_ptr as well).
*/
......@@ -524,14 +523,14 @@ namespace Dune
}
//! obtain constant reference to the coordinate function
GridFunction const& gridFunction () const
{
return *gridFunction_;
GridFunction const& gridFunction () const
{
return *gridFunction_;
}
GridFunction& gridFunction ()
{
return *gridFunction_;
GridFunction& gridFunction ()
{
return *gridFunction_;
}
bool useGeometryCaching () const
......@@ -585,7 +584,7 @@ namespace Dune
//! Deduction guide for CurvedSurfaceGrid from a grid-functions or callables
template< class HostGrid, class GF, int ORDER = -1 >
CurvedSurfaceGrid (HostGrid& hostGrid, GF&& gridFunction, std::integral_constant<int,ORDER> = {})
CurvedSurfaceGrid (HostGrid& hostGrid, GF&& gridFunction, std::integral_constant<int,ORDER> = {})
-> CurvedSurfaceGrid<GridFunctionOf_t<HostGrid,std::decay_t<GF>>,ORDER>;
......
......@@ -19,10 +19,6 @@
#include <dune/curvedsurfacegrid/localgeometrywrapper.hh>
#include <dune/curvedsurfacegrid/gridfunctions/gridfunction.hh>
#include <dune/functions/common/functionconcepts.hh>
#include <dune/functions/common/signature.hh>
#include <dune/functions/gridfunctions/localderivativetraits.hh>
namespace Dune
{
......@@ -33,22 +29,11 @@ namespace Dune
{
template< class GF >
struct DimRange
{
using SigTraits = Functions::SignatureTraits<GF>;
using Range = typename SigTraits::RawRange;
static const int value = Range::size();
};
template< class GF, class LF >
struct DifferentiableLocalFunction
{
using EntitySet = typename GF::EntitySet;
using LocalContext = typename EntitySet::Element;
using Range = typename Functions::SignatureTraits<GF>::Range;
using LocalSignature = Range(typename EntitySet::LocalCoordinate);
static const bool value = Functions::Concept::isDifferentiableLocalFunction<LF,LocalSignature,LocalContext,Functions::LocalDerivativeTraits<EntitySet>::template Traits>();
using Range = std::result_of_t<GF(typename EntitySet::GlobalCoordinate)>;
using RawRange = std::decay_t<Range>;
static const int value = RawRange::size();
};
// GridFamily
......@@ -65,9 +50,6 @@ namespace Dune
using GridFunction = GF;
using LocalFunction = std::decay_t<decltype(localFunction(std::declval<GF const&>()))>;
static const bool differentiableLocalFunction = DifferentiableLocalFunction<GF, LocalFunction>::value;
static_assert((differentiableLocalFunction || order>0), "Either provide a differentiable GridFunction or set ORDER > 0");
using ctype = typename HostGrid::ctype;
static const int dimension = HostGrid::dimension;
......
install(FILES
install(FILES
analyticgridfunction.hh
discretegridviewfunction.hh
ellipsoidgridfunction.hh
gridentityset.hh
gridfunction.hh
normalgridviewfunction.hh
......
......@@ -3,11 +3,11 @@
#ifndef DUNE_CURVED_SURFACE_GRID_ANALYTIC_GRIDFUNCTION_HH
#define DUNE_CURVED_SURFACE_GRID_ANALYTIC_GRIDFUNCTION_HH
#include <optional>
#include <type_traits>
#include <utility>
#include <dune/common/typeutilities.hh>
#include <dune/common/std/optional.hh>
#include "gridentityset.hh"
......@@ -102,8 +102,8 @@ namespace Dune
Functor f_;
// some caches
Std::optional<LocalContext> localContext_;
Std::optional<Geometry> geometry_;
std::optional<LocalContext> localContext_;
std::optional<Geometry> geometry_;
};
//! Derivative of a \ref LocalAnalyticGridFunction
......
......@@ -6,6 +6,10 @@
#include <array>
#include <vector>
#if !HAVE_DUNE_FUNCTIONS
#error "Need dune-functions for the definition of the DiscreteGridViewFunction"
#endif
#include <dune/common/fvector.hh>
#include <dune/functions/backends/istlvectorbackend.hh>
#include <dune/functions/common/defaultderivativetraits.hh>
......@@ -54,7 +58,7 @@ namespace Dune
/**
* \param gridView Function that can be evaluated at global coordinates of the \ref Grid
* \param order Polynomial order of the local parametrization [ORDER]
*
*
* \tparam ORDER Polynomial order of the local parametrization [-1]
* \tparam T Range-type of the basis and coefficient value-type [double]
**/
......@@ -152,7 +156,7 @@ namespace Dune
//! Evaluate coordinates in local coordinates
//! by interpolation of stored coords in \ref localCoords_.
Range operator() (const Domain& local) const
Range operator() (const Domain& local) const
{
static_assert(derivativeOrder < 2, "Higher-order derivatives not implemented");
......@@ -170,7 +174,7 @@ namespace Dune
}
private:
DerivativeRange<0> evaluateFunction (const Domain& local) const
DerivativeRange<0> evaluateFunction (const Domain& local) const
{
assert(bound_);
......@@ -188,7 +192,7 @@ namespace Dune
return x;
}
DerivativeRange<1> evaluateJacobian (const Domain& local) const
DerivativeRange<1> evaluateJacobian (const Domain& local) const
{
assert(bound_);
......@@ -271,12 +275,12 @@ namespace Dune
return entitySet_;
}
const Basis& basis () const
const Basis& basis () const
{
return basis_;
}
const VectorType& coefficients () const
const VectorType& coefficients () const
{
return coords_;
}
......
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_CURVED_SURFACE_GRID_SPHERE_GRIDFUNCTION_HH
#define DUNE_CURVED_SURFACE_GRID_SPHERE_GRIDFUNCTION_HH
#ifndef DUNE_CURVED_SURFACE_GRID_ELLIPSOID_GRIDFUNCTION_HH
#define DUNE_CURVED_SURFACE_GRID_ELLIPSOID_GRIDFUNCTION_HH
#include <type_traits>
#include <cmath>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/math.hh>
#include <dune/functions/common/defaultderivativetraits.hh>
#include "analyticgridfunction.hh"
namespace Dune
......@@ -16,6 +16,9 @@ namespace Dune
template< class T >
class EllipsoidProjection
{
using Domain = FieldVector<T,3>;
using Jacobian = FieldMatrix<T,3,3>;
T a_;
T b_;
T c_;
......@@ -29,58 +32,35 @@ namespace Dune
{}
//! project the coordinate to the ellipsoid
// NOTE: This is not a closes-point projection, but a spherical-coordinate projection
template< class Domain >
// NOTE: This is not a closest-point projection, but a spherical-coordinate projection
Domain operator() (const Domain& X) const
{
using std::sin; using std::cos;
auto [phi,theta] = angles(X);
return {a_*cos(phi)*sin(theta), b_*sin(phi)*sin(theta), c_*cos(theta)};
using std::sqrt;
T x = X[0], y = X[1], z = X[2];
T nrm = sqrt(x*x + y*y + z*z);
return { a_*x/nrm , b_*y/nrm , c_*z/nrm };
}
//! derivative of the projection
friend auto derivative (const EllipsoidProjection& ellipsoid)
{
return [a=ellipsoid.a_,b=ellipsoid.b_,c=ellipsoid.c_](auto const& X)
return [a=ellipsoid.a_,b=ellipsoid.b_,c=ellipsoid.c_](Domain const& X)
{
using std::sqrt;
using Domain = std::decay_t<decltype(X)>;
using DerivativeTraits = Functions::DefaultDerivativeTraits<Domain(Domain)>;
typename DerivativeTraits::Range out;
T x = X[0], y = X[1], z = X[2];
T x2 = x*x, y2 = y*y, z2 = z*z;
T x5 = x2*x2*x;
T nrm0 = x2 + y2;
T nrm1 = x2 + y2 + z2;
T nrm2 = sqrt(nrm0/nrm1);
T nrm3 = sqrt(nrm0/x2);
T nrm4 = sqrt(nrm1)*nrm1;
T nrm5 = sqrt(nrm0)*nrm4;
return {
{
a*x*nrm3*(y2 + z2)/nrm5 ,
-b*y*nrm0/(nrm3*nrm5) ,
-c*z*x/nrm4
},
{
-a*x2*y*nrm3/nrm5 ,
b*nrm0*nrm0*nrm0*(x2 + z2)/(x5*power(nrm3, 5)*nrm5) ,
-c*y*z/nrm4
},
{
-a*z*nrm0/(nrm3*nrm5) ,
-b*y*z*nrm0/(x*nrm3*nrm5) ,
c*(x2 + y2)/nrm4
}
T nrm2 = x2 + y2 + z2;
T nrm3 = sqrt(nrm2)*nrm2;
return Jacobian{
{ a*(y2 + z2)/nrm3 , -b*y*x/nrm3 , -c*z*x/nrm3 },
{-a*x*y/nrm3 , b*(x2 + z2)/nrm3 , -c*z*y/nrm3 },
{-a*x*z/nrm3 , -b*y*z/nrm3 , c*(x2 + y2)/nrm3 }
};
};
}
//! Normal vector
template< class Domain >
Domain normal (const Domain& X) const
{
using std::sqrt;
......@@ -88,39 +68,39 @@ namespace Dune
T a2 = a_*a_, b2 = b_*b_, c2 = c_*c_;
auto div = sqrt(b2*b2*c2*c2*x*x + a2*a2*c2*c2*y*y + a2*a2*b2*b2*z*z);
return {b2*c2*x/div, a2*c2*y/div, a2*b2*z/div};
return { b2*c2*x/div , a2*c2*y/div , a2*b2*z/div };
}
//! Mean curvature
template< class Domain >
T mean_curvature (const Domain& X) const
{
using std::sqrt; using std::abs;
T x = X[0], y = X[1], z = X[2];
T x2 = x*x, y2 = y*y, z2 = z*z;
T a2 = a_*a_, b2 = b_*b_, c2 = c_*c_;
auto div = 2*a2*b2*c2*power(sqrt(x*x/(a2*a2) + y*y/(b2*b2) + z*z/(c2*c2)), 3);
return abs(x*x + y*y + z*z - a2 - b2 - c2)/div;
auto div = 2*a2*b2*c2 * power(sqrt(x2/(a2*a2) + y2/(b2*b2) + z2/(c2*c2)), 3);
return abs(x2 + y2 + z2 - a2 - b2 - c2)/div;
}
//! Gaussian curvature
template< class Domain >
T gauss_curvature (const Domain& X) const
T gauss_curvature (const Domain& X) const
{
T x = X[0], y = X[1], z = X[2];
T x2 = x*x, y2 = y*y, z2 = z*z;
T a2 = a_*a_, b2 = b_*b_, c2 = c_*c_;
auto div = a2*b2*c2*power(x*x/(a2*a2) + y*y/(b2*b2) + z*z/(c2*c2), 2);
auto div = a2*b2*c2 * power(x2/(a2*a2) + y2/(b2*b2) + z2/(c2*c2), 2);
return T(1)/div;
}
private:
FieldVector<T,2> angles (Domain x) const
std::array<T,2> angles (Domain x) const
{
using std::acos; using std::atan2;
x /= x.two_norm();
return {atan2(x[1], x[0]), acos(x[2])};
return { atan2(x[1], x[0]), acos(x[2]) };
}
};
......@@ -133,4 +113,4 @@ namespace Dune
} // end namespace Dune
#endif // DUNE_CURVED_SURFACE_GRID_SPHERE_GRIDFUNCTION_HH
#endif // DUNE_CURVED_SURFACE_GRID_ELLIPSOID_GRIDFUNCTION_HH
......@@ -4,6 +4,7 @@
#define DUNE_CURVED_SURFACE_GRID_GRIDFUNCTION_HH
#include <type_traits>
#include <dune/curvedsurfacegrid/concepts.hh>
#include <dune/curvedsurfacegrid/gridfunctions/analyticgridfunction.hh>
......@@ -32,6 +33,7 @@ namespace Dune
template< class GF >
using GridOf_t = typename Impl::GridOf<std::decay_t<decltype(std::declval<GF>().entitySet())>>::type;