Skip to content

Commit 6880440

Browse files
kendrickshepherdkendrickshepherd
committed
Made all functionality go from parent domain to spatial domain rather than parametric domain to spatial domain for the piola transformations
1 parent b3a62f2 commit 6880440

2 files changed

Lines changed: 98 additions & 22 deletions

File tree

‎src/eval/SplineSpaceEvaluator.cpp‎

Lines changed: 25 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -129,17 +129,14 @@ namespace eval
129129
return jac.determinant();
130130
}
131131

132-
Eigen::VectorXd paramToSpatialGradDeterminant( const SplineSpaceEvaluator& geom_evals, const Eigen::MatrixXd& cpts )
132+
Eigen::VectorXd gradDeterminant( const Eigen::MatrixXd& J, const Eigen::MatrixXd& H, const size_t param_dim, const size_t spatial_dim )
133133
{
134-
const Eigen::MatrixXd J = geom_evals.evaluateParamToSpatialJacobian( cpts );
135-
const Eigen::MatrixXd H = geom_evals.evaluateParamToSpatialHessian( cpts );
136-
137-
if( geom_evals.splineSpace().basisComplex().parametricAtlas().cmap().dim() == 2 and cpts.rows() == 2 )
134+
if( param_dim == 2 and spatial_dim == 2 )
138135
{
139136
return Eigen::Vector2d( H( 0, 0 ) * J( 1, 1 ) + J( 0, 0 ) * H( 1, 1 ) - H( 1, 0 ) * J( 0, 1 ) - J( 1, 0 ) * H( 0, 1 ),
140137
H( 0, 1 ) * J( 1, 1 ) + J( 0, 0 ) * H( 1, 2 ) - H( 1, 1 ) * J( 0, 1 ) - J( 1, 0 ) * H( 0, 2 ) );
141138
}
142-
else if( geom_evals.splineSpace().basisComplex().parametricAtlas().cmap().dim() == 3 )
139+
else if( param_dim == 3 )
143140
{
144141
const double& x_s = J(0,0), x_t = J(0,1), x_u = J(0,2);
145142
const double& y_s = J(1,0), y_t = J(1,1), y_u = J(1,2);
@@ -179,6 +176,14 @@ namespace eval
179176
}
180177
}
181178

179+
Eigen::VectorXd parentToSpatialGradDeterminant( const SplineSpaceEvaluator& geom_evals, const Eigen::MatrixXd& cpts )
180+
{
181+
return gradDeterminant( geom_evals.evaluateJacobian( cpts ),
182+
geom_evals.evaluateHessian( cpts ),
183+
geom_evals.splineSpace().basisComplex().parametricAtlas().cmap().dim(),
184+
cpts.rows() );
185+
}
186+
182187
Eigen::MatrixXd piolaTransformedH1FirstDerivatives( const SplineSpaceEvaluator& scalar_evals,
183188
const SplineSpaceEvaluator& geom_evals,
184189
const Eigen::MatrixXd& cpts )
@@ -192,7 +197,7 @@ namespace eval
192197
const SplineSpaceEvaluator& geom_evals,
193198
const Eigen::MatrixXd& cpts )
194199
{
195-
const Eigen::MatrixXd jac_inverse = geom_evals.evaluateParamToSpatialJacobian( cpts ).inverse();
200+
const Eigen::MatrixXd jac_inverse = geom_evals.evaluateJacobian( cpts ).inverse();
196201
return vec_evals.evaluateBasis() * jac_inverse;
197202
}
198203

@@ -210,11 +215,11 @@ namespace eval
210215
const size_t param_dim = geom_evals.splineSpace().basisComplex().parametricAtlas().cmap().dim();
211216
const size_t spatial_dim = cpts.rows();
212217
const size_t n_funcs = vector_basis.cols();
213-
const Eigen::MatrixXd jac = geom_evals.evaluateParamToSpatialJacobian( cpts );
218+
const Eigen::MatrixXd jac = geom_evals.evaluateJacobian( cpts );
214219
const auto inverse_transpose_jacobian = jac.inverse().transpose();
215220

216221
const auto modified_hessian = [&geom_evals, &cpts, &param_dim, &spatial_dim]() {
217-
const Eigen::MatrixXd hess = geom_evals.evaluateParamToSpatialHessian( cpts );
222+
const Eigen::MatrixXd hess = geom_evals.evaluateHessian( cpts );
218223
Eigen::MatrixXd out( spatial_dim * param_dim, param_dim );
219224
if( param_dim == 2 )
220225
// FIXME: Make this work for 3d spatial/2d manifold
@@ -249,8 +254,7 @@ namespace eval
249254
}();
250255

251256
return ( ( inverse_transpose_jacobian *
252-
vec_evals.evaluateFirstDerivativesFromParamToSpatial().transpose().reshaped( param_dim,
253-
n_funcs * param_dim ) )
257+
vec_evals.evaluateFirstDerivatives().transpose().reshaped( param_dim, n_funcs * param_dim ) )
254258
.reshaped( spatial_dim * param_dim, n_funcs ) +
255259
second_term )
256260
.transpose() *
@@ -261,7 +265,7 @@ namespace eval
261265
const SplineSpaceEvaluator& geom_evals,
262266
const Eigen::MatrixXd& cpts )
263267
{
264-
const auto jac = geom_evals.evaluateParamToSpatialJacobian( cpts );
268+
const auto jac = geom_evals.evaluateJacobian( cpts );
265269
const double det = determinant( jac );
266270
return 1.0 / det * vec_evals.evaluateBasis() * jac.transpose();
267271
}
@@ -270,7 +274,7 @@ namespace eval
270274
const SplineSpaceEvaluator& geom_evals,
271275
const Eigen::MatrixXd& cpts )
272276
{
273-
const Eigen::MatrixXd jac = geom_evals.evaluateParamToSpatialJacobian( cpts );
277+
const Eigen::MatrixXd jac = geom_evals.evaluateJacobian( cpts );
274278
const double det_inverse = 1.0 / determinant( jac );
275279
const Eigen::MatrixXd vector_basis = vec_evals.evaluateBasis().transpose();
276280
const size_t n_funcs = vector_basis.cols();
@@ -282,11 +286,11 @@ namespace eval
282286

283287
const Eigen::MatrixXd first_term =
284288
-det_inverse * det_inverse *
285-
( ( jac * vector_basis ).reshaped() * paramToSpatialGradDeterminant( geom_evals, cpts ).transpose() )
289+
( ( jac * vector_basis ).reshaped() * parentToSpatialGradDeterminant( geom_evals, cpts ).transpose() )
286290
.reshaped( spatial_dim * param_dim, n_funcs );
287291

288292
const auto modified_hessian = [&geom_evals, &cpts, &param_dim, &spatial_dim]() {
289-
const Eigen::MatrixXd hess = geom_evals.evaluateParamToSpatialHessian( cpts );
293+
const Eigen::MatrixXd hess = geom_evals.evaluateHessian( cpts );
290294
Eigen::MatrixXd out( spatial_dim * param_dim, param_dim );
291295
if( param_dim == 2 )
292296
// FIXME: Make this work for 3d spatial/2d manifold
@@ -312,7 +316,7 @@ namespace eval
312316
const Eigen::MatrixXd second_term = det_inverse * modified_hessian() * vector_basis;
313317

314318
const Eigen::MatrixXd third_term =
315-
( det_inverse * jac * vec_evals.evaluateFirstDerivativesFromParamToSpatial().transpose().reshaped( param_dim, n_funcs * param_dim ) )
319+
( det_inverse * jac * vec_evals.evaluateFirstDerivatives().transpose().reshaped( param_dim, n_funcs * param_dim ) )
316320
.reshaped( spatial_dim * param_dim, n_funcs );
317321

318322
return ( first_term + second_term + third_term ).transpose() *
@@ -323,17 +327,17 @@ namespace eval
323327
const SplineSpaceEvaluator& geom_evals,
324328
const Eigen::MatrixXd& cpts )
325329
{
326-
return 1.0 / determinant( geom_evals.evaluateParamToSpatialJacobian( cpts ) ) * bivec_evals.evaluateBasis();
330+
return 1.0 / determinant( geom_evals.evaluateJacobian( cpts ) ) * bivec_evals.evaluateBasis();
327331
}
328332
Eigen::MatrixXd piolaTransformedL2FirstDerivatives( const SplineSpaceEvaluator& bivec_evals,
329333
const SplineSpaceEvaluator& geom_evals,
330334
const Eigen::MatrixXd& cpts )
331335
{
332-
const auto jac = geom_evals.evaluateParamToSpatialJacobian( cpts );
336+
const auto jac = geom_evals.evaluateJacobian( cpts );
333337
const double det_inverse = 1.0 / determinant( jac );
334338
const auto jac_inverse_transpose = jac.inverse().transpose();
335-
return ( det_inverse * jac_inverse_transpose * bivec_evals.evaluateFirstDerivativesFromParamToSpatial().transpose() -
336-
det_inverse * det_inverse * jac_inverse_transpose * paramToSpatialGradDeterminant( geom_evals, cpts ) *
339+
return ( det_inverse * jac_inverse_transpose * bivec_evals.evaluateFirstDerivatives().transpose() -
340+
det_inverse * det_inverse * jac_inverse_transpose * parentToSpatialGradDeterminant( geom_evals, cpts ) *
337341
bivec_evals.evaluateBasis().transpose() )
338342
.transpose() * jac.inverse(); // transform from ds denominator to dx
339343
}
@@ -470,4 +474,4 @@ namespace eval
470474
{
471475
throw std::runtime_error( "evaluateFirstDerivativesFromParamToSpatial not implemented for NURBSSpaceEvaluator" );
472476
}
473-
}
477+
}

‎test/SplineSpaceEvaluator_test.cpp‎

Lines changed: 73 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -95,6 +95,78 @@ TEST_CASE( "TP Spline space evaluation" )
9595
} );
9696
}
9797

98+
TEST_CASE( "Piola transforms are invariant to parent-to-parametric interval scaling" )
99+
{
100+
const double ptol = 1e-10;
101+
const size_t degree = 1;
102+
103+
const auto build_discretization = [&]( const KnotVector& kv_s, const KnotVector& kv_t ) {
104+
const TPSplineSpace h1 = buildBSpline( { kv_s, kv_t }, { degree, degree } );
105+
const auto hdiv_bc = std::make_shared<const VectorConformingBasisComplex>( h1.basisComplexPtr() );
106+
const VectorConformingTPSplineSpace hdiv( hdiv_bc, h1 );
107+
const auto l2_bc =
108+
std::make_shared<const TPBasisComplex>( h1.basisComplex().parametricAtlasPtr(),
109+
hdiv.reducedDegree1dBases().at( 0 )->basisComplexPtr(),
110+
hdiv.reducedDegree1dBases().at( 1 )->basisComplexPtr() );
111+
const TPSplineSpace l2( l2_bc, hdiv.reducedDegree1dBases().at( 0 ), hdiv.reducedDegree1dBases().at( 1 ) );
112+
return std::tuple<TPSplineSpace, VectorConformingTPSplineSpace, TPSplineSpace>( h1, hdiv, l2 );
113+
};
114+
115+
const KnotVector kv_unit( { 0, 0, 1, 1 }, ptol );
116+
const KnotVector kv_double( { 0, 0, 2, 2 }, ptol );
117+
118+
auto [h1_unit, hdiv_unit, l2_unit] = build_discretization( kv_unit, kv_unit );
119+
auto [h1_scaled, hdiv_scaled, l2_scaled] = build_discretization( kv_double, kv_unit );
120+
121+
eval::SplineSpaceEvaluator h1_eval_unit( h1_unit, 2 );
122+
eval::SplineSpaceEvaluator hdiv_eval_unit( hdiv_unit, 1 );
123+
eval::SplineSpaceEvaluator l2_eval_unit( l2_unit, 1 );
124+
125+
eval::SplineSpaceEvaluator h1_eval_scaled( h1_scaled, 2 );
126+
eval::SplineSpaceEvaluator hdiv_eval_scaled( hdiv_scaled, 1 );
127+
eval::SplineSpaceEvaluator l2_eval_scaled( l2_scaled, 1 );
128+
129+
const topology::Face elem( topology::Dart( 0 ) );
130+
const param::ParentPoint pt( param::cubeDomain( 2 ), Eigen::Vector2d( 0.37, 0.61 ), { false, false, false, false } );
131+
132+
h1_eval_unit.localizeElement( elem );
133+
hdiv_eval_unit.localizeElement( elem );
134+
l2_eval_unit.localizeElement( elem );
135+
h1_eval_scaled.localizeElement( elem );
136+
hdiv_eval_scaled.localizeElement( elem );
137+
l2_eval_scaled.localizeElement( elem );
138+
139+
h1_eval_unit.localizePoint( pt );
140+
hdiv_eval_unit.localizePoint( pt );
141+
l2_eval_unit.localizePoint( pt );
142+
h1_eval_scaled.localizePoint( pt );
143+
hdiv_eval_scaled.localizePoint( pt );
144+
l2_eval_scaled.localizePoint( pt );
145+
146+
Eigen::MatrixX2d unit_geom( 4, 2 );
147+
unit_geom << 0.0, 0.0,
148+
1.0, 0.0,
149+
0.0, 1.0,
150+
1.0, 1.0;
151+
const Eigen::MatrixX2d scaled_geom = unit_geom;
152+
153+
CHECK( ( eval::piolaTransformedH1FirstDerivatives( h1_eval_unit, h1_eval_unit, unit_geom.transpose() ) -
154+
eval::piolaTransformedH1FirstDerivatives( h1_eval_scaled, h1_eval_scaled, scaled_geom.transpose() ) )
155+
.norm() < 1e-12 );
156+
CHECK( ( eval::piolaTransformedHDivBasis( hdiv_eval_unit, h1_eval_unit, unit_geom.transpose() ) -
157+
eval::piolaTransformedHDivBasis( hdiv_eval_scaled, h1_eval_scaled, scaled_geom.transpose() ) )
158+
.norm() < 1e-12 );
159+
CHECK( ( eval::piolaTransformedHDivFirstDerivatives( hdiv_eval_unit, h1_eval_unit, unit_geom.transpose() ) -
160+
eval::piolaTransformedHDivFirstDerivatives( hdiv_eval_scaled, h1_eval_scaled, scaled_geom.transpose() ) )
161+
.norm() < 1e-12 );
162+
CHECK( ( eval::piolaTransformedL2Basis( l2_eval_unit, h1_eval_unit, unit_geom.transpose() ) -
163+
eval::piolaTransformedL2Basis( l2_eval_scaled, h1_eval_scaled, scaled_geom.transpose() ) )
164+
.norm() < 1e-12 );
165+
CHECK( ( eval::piolaTransformedL2FirstDerivatives( l2_eval_unit, h1_eval_unit, unit_geom.transpose() ) -
166+
eval::piolaTransformedL2FirstDerivatives( l2_eval_scaled, h1_eval_scaled, scaled_geom.transpose() ) )
167+
.norm() < 1e-12 );
168+
}
169+
98170
Eigen::MatrixXd quadballControlPoints()
99171
{
100172
// See table 4 of "Tiling the Sphere with Rational Bezier Patches" (James E. Cobb, 1988)
@@ -283,4 +355,4 @@ TEST_CASE( "NURBS Spline space evaluation" )
283355

284356
return true;
285357
} );
286-
}
358+
}

0 commit comments

Comments
 (0)