From 484f8724e9afc268a10c72ec3bdd4ead4a39b1f9 Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
Date: Wed, 4 Mar 2020 15:18:59 +0100
Subject: [PATCH] add corner implementation
---
dune/curvedgeometry/curvedgeometry.hh | 57 +++++++++++--------
.../curvedgeometrydatacollector.hh | 5 +-
2 files changed, 37 insertions(+), 25 deletions(-)
diff --git a/dune/curvedgeometry/curvedgeometry.hh b/dune/curvedgeometry/curvedgeometry.hh
index ce397fa..4bdfc7b 100644
--- a/dune/curvedgeometry/curvedgeometry.hh
+++ b/dune/curvedgeometry/curvedgeometry.hh
@@ -59,7 +59,7 @@ namespace Dune
* This structure provides the default values.
*
* \tparam ct coordinate type
- * \tparam LFECache A LocalFiniteElementVariantCache implementation,
+ * \tparam LFECache A LocalFiniteElementVariantCache implementation,
* e.g. LagrangeLocalFiniteElementCache
*/
template <class ct, class LFECache>
@@ -107,8 +107,10 @@ namespace Dune
/// type of local coordinates
using LocalCoordinate = FieldVector<ctype, mydimension>;
+
/// type of global coordinates
using GlobalCoordinate = FieldVector<ctype, coorddimension>;
+
/// type of volume
using Volume = ctype;
@@ -118,12 +120,9 @@ namespace Dune
/// type of jacobian inverse transposed
using JacobianInverseTransposed = FieldMatrix<ctype, coorddimension, mydimension>;
- protected:
- using ReferenceElements = Dune::ReferenceElements<ctype, mydimension>;
-
public:
/// type of reference element
- using ReferenceElement = typename ReferenceElements::ReferenceElement;
+ using ReferenceElement = typename ReferenceElements<ctype, mydimension>::ReferenceElement;
/// Parametrization of the geometry
using Traits = TraitsType;
@@ -137,6 +136,7 @@ namespace Dune
using LocalBasisTraits = typename LocalBasis::Traits;
protected:
+ // Internal constructor. Must not be called directly since more initialization is necessary
CurvedGeometry (const ReferenceElement& refElement)
: refElement_(refElement)
, localFECache_()
@@ -150,9 +150,8 @@ namespace Dune
* \param[in] refElement reference element for the geometry
* \param[in] vertices vertices to store internally
*
- * \note The type of vertices is actually a template argument.
- * It is only required that the internal vertex storage can be
- * constructed from this object.
+ * \note The vertices are stored internally, so if possible move an external vertex storedge
+ * to this constructor
*/
CurvedGeometry (const ReferenceElement& refElement, std::vector<GlobalCoordinate> vertices)
: CurvedGeometry(refElement)
@@ -233,13 +232,13 @@ namespace Dune
/// \brief obtain coordinates of the i-th corner
/**
- * \note The corners are part of the stored vertices, always stored first before
- * higher order vertices.
+ * \todo In some local-bases, the corners are part of the stored vertices, and sometimes
+ * even stored first before higher order vertices. Maybe this can be utilized somehow.
*/
GlobalCoordinate corner (int i) const
{
assert( (i >= 0) && (i < corners()) );
- return vertices_[i];
+ return global(refElement_.position(i, mydimension));
}
/// \brief obtain the centroid of the mapping's image
@@ -255,7 +254,7 @@ namespace Dune
*
* \f[ global = \sum_i v_i \psi_i(local) \f]
*
- * \param[in] local local coordinate to map
+ * \param[in] local local coordinate to map
*
* \returns corresponding global coordinate
*/
@@ -278,7 +277,7 @@ namespace Dune
*
* \return corresponding local coordinate
*
- * \note For given global coordinate y the returned local coordinate x that minimizes
+ * \note For given global coordinate `y` the returned local coordinate `x` that minimizes
* the following function over the local coordinate space spanned by the reference element.
* \code
* (global( x ) - y).two_norm()
@@ -297,6 +296,7 @@ namespace Dune
const GlobalCoordinate dglobal = global(x) - globalCoord;
const bool invertible = MatrixHelper::xTRightInvA(jacobianTransposed(x), dglobal, dx);
+ // break if jacobian is not invertible
if (!invertible)
break;
@@ -307,6 +307,7 @@ namespace Dune
if (affineMapping)
break;
+ // break if tolerance is reached.
if (dx.two_norm2() < tolerance)
break;
}
@@ -315,7 +316,8 @@ namespace Dune
return x;
}
- /// \brief Construct a global coordinate normal of the curvilinear element evaluated at a given local coordinate
+ /// \brief Construct a global coordinate normal of the curvilinear element evaluated at
+ /// a given local coordinate
/**
* \note Implemented for codim=1 entities only, i.e. edges in 2D and faces in 3D
**/
@@ -338,9 +340,8 @@ namespace Dune
*
* \returns the integration element \f$\mu(x)\f$.
*
- * \note For affine mappings, it is more efficient to call
- * jacobianInverseTransposed before integrationElement, if both
- * are required.
+ * \todo Implement caching of the jacobianTransposed matrix, since it is used fir the
+ * integrationElement and the jacobianInverseTransposed
*/
ctype integrationElement (const LocalCoordinate& local) const
{
@@ -350,10 +351,12 @@ namespace Dune
/// \brief Obtain the volume of the mapping's image
/**
* Calculates the volume of the entity by numerical integration
+ *
+ * \todo Check the quadrature order of the volume integral.
*/
Volume volume () const
{
- const int p = 2*localBasis().order(); // TODO: needs to be checked
+ const int p = 2*localBasis().order();
const auto& quadRule = QuadratureRules<ctype, mydimension>::rule(type(), p);
Volume vol(0);
for (auto const& qp : quadRule)
@@ -366,10 +369,7 @@ namespace Dune
/**
* \param[in] local local coordinate to evaluate Jacobian in
*
- * \returns a reference to the transposed of the Jacobian
- *
- * \note The returned reference is reused on the next call to
- * JacobianTransposed, destroying the previous value.
+ * \returns the matrix corresponding to the transposed of the Jacobian
*/
JacobianTransposed jacobianTransposed (const LocalCoordinate& local) const
{
@@ -410,16 +410,19 @@ namespace Dune
}
protected:
- ReferenceElement refElement () const
+ // the internal stored reference element
+ const ReferenceElement& refElement () const
{
return refElement_;
}
+ // the local basis of the stored local finite-element
const LocalBasis& localBasis() const
{
return localFE_.localBasis();
}
+ // normal vector to an edge line-element
GlobalCoordinate normalEdge (const LocalCoordinate& local, const JacobianTransposed& J) const
{
GlobalCoordinate res{
@@ -429,6 +432,7 @@ namespace Dune
return res /= res.two_norm();
}
+ // normal vector to a triangle face element
GlobalCoordinate normalTriangle (const LocalCoordinate& local, const JacobianTransposed& J) const
{
GlobalCoordinate res{
@@ -440,13 +444,20 @@ namespace Dune
}
private:
+ //! Reference of the geometry
ReferenceElement refElement_;
+
+ //! Cache to construct finite elements
LocalFECache localFECache_;
+
+ //! A cached local finite-element
LocalFiniteElement const& localFE_;
+ //! The (lagrange) coefficients of the interpolating geometry
std::vector<GlobalCoordinate> vertices_;
};
+
/// \brief curved geometry parametrized with lagrange basis functions
template <class ctype, int mydim, int cdim, int order>
using LagrangeCurvedGeometry = CurvedGeometry<ctype,mydim,cdim,
diff --git a/dune/curvedgeometry/curvedgeometrydatacollector.hh b/dune/curvedgeometry/curvedgeometrydatacollector.hh
index 67b7c4a..ad90e1e 100644
--- a/dune/curvedgeometry/curvedgeometrydatacollector.hh
+++ b/dune/curvedgeometry/curvedgeometrydatacollector.hh
@@ -18,10 +18,11 @@ namespace Dune
// and write nonlinear vtk elements
template <class GridView, class Geometry>
class CurvedGeometryDataCollector
- : public UnstructuredDataCollectorInterface<GridView, CurvedGeometryDataCollector<GridView,Geometry>, Partitions::All>
+ : public UnstructuredDataCollectorInterface<GridView,
+ CurvedGeometryDataCollector<GridView,Geometry>, Partitions::All>
{
using Self = CurvedGeometryDataCollector;
- using Super = UnstructuredDataCollectorInterface<GridView, CurvedGeometryDataCollector<GridView,Geometry>, Partitions::All>;
+ using Super = UnstructuredDataCollectorInterface<GridView, Self, Partitions::All>;
using LocalCoord = typename GridView::template Codim<0>::Entity::Geometry::LocalCoordinate;
using GlobalCoord = typename GridView::template Codim<0>::Entity::Geometry::GlobalCoordinate;
--
GitLab