Getting the Parametric Points

Author: Raminpahnabi

src/api/SplinesAPI.cpp

@@ -161,6 +161,11 @@ PYBIND11_MODULE( splines, m )
             []( const api::NavierStokesDiscretization& nsd ) { return nsd.getH1().evaluateManifold( nsd.controlPoints() ); },
             "Evaluates the spatial position of the spline geometry at the parametric position from the latest calls to "
             "localizeElement and localizePoint." )
+        .def(
+            "parametricPoint",
+            []( const api::NavierStokesDiscretization& nsd ) { return nsd.getH1().evaluateParametricPoint(); },
+            "Evaluates the parametric coordinates (knot-space) of the current quadrature point from the latest calls to "
+            "localizeElement and localizePoint." )
         .def(
             "jacobian",
             []( const api::NavierStokesDiscretization& nsd ) { return nsd.getH1().evaluateJacobian( nsd.controlPoints() ); },

src/eval/SplineSpaceEvaluator.cpp

@@ -18,12 +18,21 @@ namespace eval
         mExOp = mSpline.extractionOperator( c );
         mCurrentCell.emplace( c );
         mParametricLengths = mSpline.basisComplex().parametricAtlas().parametricLengths( c );
+        mParametricStart = mSpline.basisComplex().parametricAtlas().parametricStarts( c );
     }
 
     void SplineSpaceEvaluator::localizePoint( const param::ParentPoint& ppt )
     {
         if( not mCurrentCell.has_value() ) throw std::runtime_error( "Must localize cell before locaizing point" );
         mLocalEvals.emplace( ParentBasisEval( mSpline.basisComplex().parentBasis( mCurrentCell.value() ), ppt, mNumDerivs ) );
+        mLocalizedParentCoords.emplace( ppt.mPoint.cast<double>() );
+    }
+
+    Eigen::VectorXd SplineSpaceEvaluator::evaluateParametricPoint() const
+    {
+        if( not mLocalizedParentCoords.has_value() ) throw std::runtime_error( "Must localize point before evaluating parametric point" );
+        const size_t param_dim = mSpline.basisComplex().parametricAtlas().cmap().dim();
+        return mParametricStart.head( param_dim ) + mLocalizedParentCoords->cwiseProduct( mParametricLengths.head( param_dim ) );
     }
 
     Eigen::VectorXd SplineSpaceEvaluator::evaluateManifold( const Eigen::MatrixXd& cpts ) const

src/eval/SplineSpaceEvaluator.hpp

@@ -31,6 +31,8 @@ namespace eval
 
         virtual Eigen::MatrixXd evaluateFirstDerivativesFromParamToSpatial() const;
 
+        virtual Eigen::VectorXd evaluateParametricPoint() const;
+
         size_t numDerivatives() const { return mNumDerivs; }
 
         const basis::SplineSpace& splineSpace() const { return mSpline; }
@@ -40,9 +42,11 @@ namespace eval
         const basis::SplineSpace& mSpline;
         std::optional<topology::Cell> mCurrentCell;
         Vector6dMax mParametricLengths;
+        Vector6dMax mParametricStart;
         std::vector<basis::FunctionId> mConnect;
         Eigen::MatrixXd mExOp;
         std::optional<ParentBasisEval> mLocalEvals;
+        std::optional<Eigen::VectorXd> mLocalizedParentCoords;
     };
 
     Eigen::MatrixXd piolaTransformedH1FirstDerivatives( const SplineSpaceEvaluator& h1_evals, const SplineSpaceEvaluator& geom_evals, const Eigen::MatrixXd& cpts );

src/param/HierarchicalTPParametricAtlas.cpp

@@ -158,4 +158,12 @@ Vector6dMax HierarchicalTPParametricAtlas::parametricLengths( const topology::Ce
         throw std::invalid_argument( "HierarchicalTPParametricAtlas::parametricLengths only works for elements." );
     const auto [ level, level_d ] = mMap->unrefinedAncestorDartOfCell( c );
     return mRefinementLevels.at( level )->parametricLengths( topology::Cell( level_d, c.dim() ) );
+}
+
+Vector6dMax HierarchicalTPParametricAtlas::parametricStarts( const topology::Cell& c ) const
+{
+    if( c.dim() != mMap->dim() )
+        throw std::invalid_argument( "HierarchicalTPParametricAtlas::parametricStarts only works for elements." );
+    const auto [ level, level_d ] = mMap->unrefinedAncestorDartOfCell( c );
+    return mRefinementLevels.at( level )->parametricStarts( topology::Cell( level_d, c.dim() ) );
 }

src/param/HierarchicalTPParametricAtlas.hpp

@@ -20,6 +20,7 @@ namespace param
         virtual const ParentDomain parentDomain( const topology::Cell& c ) const override;
         virtual ParentPoint parentPoint( const topology::Vertex& v ) const override;
         virtual Vector6dMax parametricLengths( const topology::Cell& c ) const override;
+        virtual Vector6dMax parametricStarts( const topology::Cell& c ) const override;
 
         const std::vector<std::shared_ptr<const TPParametricAtlas>>& refinementLevels() const { return mRefinementLevels; }
 

src/param/ParametricAtlas.hpp

@@ -11,6 +11,10 @@ namespace param
         virtual const ParentDomain parentDomain( const topology::Cell& c ) const = 0;
         virtual ParentPoint parentPoint( const topology::Vertex& v ) const = 0;
         virtual Vector6dMax parametricLengths( const topology::Cell& c ) const = 0;
+        virtual Vector6dMax parametricStarts( [[maybe_unused]] const topology::Cell& c ) const
+        {
+            throw std::runtime_error( "parametricStarts not implemented for this atlas type" );
+        }
     };
 
     BaryCoordIsZeroVec parentDomainBoundary( const ParametricAtlas& atlas, const topology::Cell& cell );

src/param/ParametricAtlas1d.cpp

@@ -29,6 +29,12 @@ namespace param
         return Eigen::Matrix<double, 1, 1>( mLengths( c.dart().id() ) );
     }
 
+    Vector6dMax ParametricAtlas1d::parametricStarts( const topology::Cell& c ) const
+    {
+        if( c.dim() != 1 ) throw std::runtime_error( "invalid cell dimension" );
+        return Eigen::Matrix<double, 1, 1>( mLengths.head( c.dart().id() ).sum() );
+    }
+
     std::vector<std::pair<topology::Cell, param::ParentPoint>> parentPointsOfParamPoints(
         const std::vector<double>& values, const param::ParametricAtlas1d& pa, const double param_tol )
     {

src/param/ParametricAtlas1d.hpp

@@ -18,6 +18,7 @@ namespace param
         virtual const ParentDomain parentDomain( const topology::Cell& c ) const override;
         virtual ParentPoint parentPoint( const topology::Vertex& v ) const override;
         virtual Vector6dMax parametricLengths( const topology::Cell& c ) const override;
+        virtual Vector6dMax parametricStarts( const topology::Cell& c ) const override;
 
         double totalLength() const { return mLengths.sum(); }
 

src/param/TPParametricAtlas.cpp

@@ -73,6 +73,14 @@ namespace param
         return ( Vector6dMax( source_lengths.size() + line_lengths.size() ) << source_lengths, line_lengths ).finished();
     }
 
+    Vector6dMax TPParametricAtlas::parametricStarts( const topology::Cell& c ) const
+    {
+        const auto [source_dart, line_dart, tp_pos] = mMap->unflatten( c.dart() );
+        const Vector6dMax source_starts = mSourceParam->parametricStarts( topology::Cell( source_dart, mSourceParam->cmap().dim() ) );
+        const Vector6dMax line_starts = mLineParam->parametricStarts( topology::Cell( line_dart, mLineParam->cmap().dim() ) );
+        return ( Vector6dMax( source_starts.size() + line_starts.size() ) << source_starts, line_starts ).finished();
+    }
+
     SmallVector<topology::Cell, 12> cornerCells( const TPParametricAtlas& atlas, const uint cell_dim )
     {
         using namespace topology;

src/param/TPParametricAtlas.hpp

@@ -19,6 +19,7 @@ namespace param
         virtual const ParentDomain parentDomain( const topology::Cell& c ) const override;
         virtual ParentPoint parentPoint( const topology::Vertex& v ) const override;
         virtual Vector6dMax parametricLengths( const topology::Cell& c ) const override;
+        virtual Vector6dMax parametricStarts( const topology::Cell& c ) const override;
 
         const ParametricAtlas& source() const { return *mSourceParam; }
         const ParametricAtlas1d& line() const { return *mLineParam; }