Skip to content

New surface point class and example using Dirichlet Trace #104

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
jdoelz wants to merge 29 commits into
base: master
Choose a base branch
from
Open

New surface point class and example using Dirichlet Trace #104

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
29 commits
Select commit Hold shift + click to select a range
43fbc25
Split storage of far field quadrature nodes from a single, large matr…
jdoelz Jun 3, 2022
8500db2
added DirichletTrace for DiscreteLinearFormOnReferenceDomain
jdoelz Jul 10, 2022
4324627
Merge remote-tracking branch 'upstream/split_ffqnodes_matrix_into_sma…
jdoelz Jan 30, 2023
0899e12
use typedef SurfacePoint consistently throughout the code
jdoelz Jan 30, 2023
68d9ab3
Typo
jdoelz Jan 30, 2023
e5b5c2d
Introduced typedef for improved compatiblity with older compilers
jdoelz Feb 1, 2023
0018349
Replaced vector representation of SurfacePoint by a wrapper class
jdoelz Feb 7, 2023
96523f4
new surface point works now, although it still seems to be a bit slow
jdoelz Feb 8, 2023
9522ab0
fix compile error
jdoelz Feb 8, 2023
6c222d9
Merge branch 'consistentsurfacepoint' into newsurfacepoint
jdoelz Feb 8, 2023
2340a3b
refactoring of shape functions, this made the code actually slower
jdoelz Feb 12, 2023
a02bf44
refactoring surfacePointUpdate
jdoelz Feb 13, 2023
95db07c
move memory allocation in updateSurfacePoint to memory allocation in …
jdoelz Feb 13, 2023
ac91125
Merge branch 'master' into DiscreteLinearFormOnReferenceDomain
jdoelz Mar 18, 2024
12de2ff
Merge remote-tracking branch 'origin/newsurfacepoint' into DiscreteLi…
jdoelz Mar 18, 2024
210bdb7
add patch number to surface point
jdoelz Mar 18, 2024
25acf64
Merge remote-tracking branch 'origin/newsurfacepoint' into DiscreteLi…
jdoelz Mar 18, 2024
2975862
debugging
jdoelz Mar 18, 2024
138fe0d
fixed bug
jdoelz Mar 19, 2024
0b0e1c9
Merge branch 'master' into newsurfacepoint
jdoelz Mar 19, 2024
62e8959
remove most of the compile errors
jdoelz Mar 19, 2024
61127df
remove minor bugs caused by merging
jdoelz Mar 20, 2024
03a26dd
Merge remote-tracking branch 'origin/newsurfacepoint_v2' into Dirichl…
jdoelz Mar 20, 2024
380c1be
style guide
jdoelz Mar 20, 2024
cdfea84
Merge remote-tracking branch 'upstream/master' into DirichletTraceOnR…
jdoelz Jun 20, 2025
ed0b2e6
Update Bembel/src/Geometry/SurfacePoint.hpp
jdoelz Jun 20, 2025
3eff3fc
fix compilation issues
jdoelz Jun 20, 2025
09281c4
Merge branch 'DirichletTraceOnReferenceDomain_v2' of github.com:temf/…
jdoelz Jun 20, 2025
f824e79
fixed style guide issue
jdoelz Jun 20, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 4 additions & 1 deletion .gitignore
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,4 +7,7 @@ eigen3
.vscode
# Important for a nice gh-pages branch
node_modules/
package-lock.json
package-lock.json
# dev
vtune
.vscode
1 change: 1 addition & 0 deletions Bembel/LinearForm
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@
#include "src/LinearForm/LinearForm.hpp"

#include "src/LinearForm/DirichletTrace.hpp"
#include "src/LinearForm/DirichletTraceOnReferenceDomain.hpp"
#include "src/LinearForm/NeumannTrace.hpp"
#include "src/LinearForm/DiscreteLinearForm.hpp"
#include "src/LinearForm/RotatedTangentialTrace.hpp"
Expand Down
1 change: 1 addition & 0 deletions Bembel/Spline
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@
#include <Eigen/Sparse>

#include "src/util/Constants.hpp"
#include "src/util/Macros.hpp"
#include "src/LinearOperator/DifferentialFormEnum.hpp"

#include "src/Spline/DeBoor.hpp"
Expand Down
2 changes: 1 addition & 1 deletion Bembel/src/AnsatzSpace/FunctionEvaluator.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -100,7 +100,7 @@ class FunctionEvaluator {

SurfacePoint sp;
ansatz_space_.get_superspace().get_geometry()[patch].updateSurfacePoint(
&sp, ref_point, 1, element.mapToReferenceElement(ref_point));
&sp, patch, ref_point, 1, element.mapToReferenceElement(ref_point));
return evaluate(element, sp);
}
Eigen::Matrix<
Expand Down
25 changes: 10 additions & 15 deletions Bembel/src/AnsatzSpace/FunctionEvaluatorEval.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,8 +29,7 @@ struct FunctionEvaluatorEval<Scalar, DifferentialForm::Continuous, LinOp> {
Scalar, Eigen::Dynamic,
getFunctionSpaceVectorDimension<DifferentialForm::Continuous>()>
&coeff) const {
auto s = p.segment<2>(0);
return coeff.transpose() * super_space.basis(s) / element.get_h();
return coeff.transpose() * super_space.basis(p.get_xi()) / element.get_h();
}
};

Expand All @@ -47,15 +46,13 @@ struct FunctionEvaluatorEval<Scalar, DifferentialForm::DivConforming, LinOp> {
Scalar, Eigen::Dynamic,
getFunctionSpaceVectorDimension<DifferentialForm::DivConforming>()>
&coeff) const {
auto s = p.segment<2>(0);
auto h = element.get_h();
auto x_f_dx = p.segment<3>(6);
auto x_f_dy = p.segment<3>(9);
double h = element.get_h();
Eigen::Matrix<typename LinearOperatorTraits<LinOp>::Scalar, Eigen::Dynamic,
1>
tangential_coefficients = coeff.transpose() * super_space.basis(s);
return (x_f_dx * tangential_coefficients(0) +
x_f_dy * tangential_coefficients(1)) /
tangential_coefficients =
coeff.transpose() * super_space.basis(p.get_xi());
return (p.get_f_dx() * tangential_coefficients(0) +
p.get_f_dy() * tangential_coefficients(1)) /
h;
}

Expand All @@ -67,14 +64,13 @@ struct FunctionEvaluatorEval<Scalar, DifferentialForm::DivConforming, LinOp> {
Scalar, Eigen::Dynamic,
getFunctionSpaceVectorDimension<DifferentialForm::DivConforming>()>
&coeff) const {
auto s = p.segment<2>(0);
auto h = element.get_h();
double h = element.get_h();
Eigen::Matrix<typename LinearOperatorTraits<LinOp>::Scalar, Eigen::Dynamic,
1>
phiPhiVec_dx = super_space.basisDx(s);
phiPhiVec_dx = super_space.basisDx(p.get_xi());
Eigen::Matrix<typename LinearOperatorTraits<LinOp>::Scalar, Eigen::Dynamic,
1>
phiPhiVec_dy = super_space.basisDy(s);
phiPhiVec_dy = super_space.basisDy(p.get_xi());
return (phiPhiVec_dx.dot(coeff.col(0)) + phiPhiVec_dy.dot(coeff.col(1))) /
h / h;
}
Expand All @@ -93,8 +89,7 @@ struct FunctionEvaluatorEval<Scalar, DifferentialForm::Discontinuous, LinOp> {
Scalar, Eigen::Dynamic,
getFunctionSpaceVectorDimension<DifferentialForm::Discontinuous>()>
&coeff) const {
auto s = p.segment<2>(0);
return coeff.transpose() * super_space.basis(s) / element.get_h();
return coeff.transpose() * super_space.basis(p.get_xi()) / element.get_h();
}
};
} // namespace Bembel
Expand Down
12 changes: 5 additions & 7 deletions Bembel/src/AnsatzSpace/Projector.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -215,18 +215,16 @@ inline _proj_info makeLocalProjectionTriplets(
// Here, we suddenly use the degree for basisevaluation, i.e.,
// maximal_polynomial_degree-1. This is confusing, but correct and tested.
{
double vals_y[Constants::MaxP + 1];
double vals_x[Constants::MaxP + 1];
for (int iy = 0; iy < maximal_polynomial_degree; ++iy) {
Bembel::Basis::ShapeFunctionHandler::evalBasis(
maximal_polynomial_degree - 1, vals_y, mask[iy]);
Eigen::VectorXd vals_y = Bembel::Basis::ShapeFunctionHandler::evalBasis(
maximal_polynomial_degree - 1, mask[iy]);
for (int ix = 0; ix < maximal_polynomial_degree; ++ix) {
Bembel::Basis::ShapeFunctionHandler::evalBasis(
maximal_polynomial_degree - 1, vals_x, mask[ix]);
Eigen::VectorXd vals_x = Bembel::Basis::ShapeFunctionHandler::evalBasis(
maximal_polynomial_degree - 1, mask[ix]);
for (int jy = 0; jy < maximal_polynomial_degree; ++jy) {
for (int jx = 0; jx < maximal_polynomial_degree; ++jx) {
system(iy * masksize + ix, jy * maximal_polynomial_degree + jx) =
vals_x[jx] * vals_y[jy];
vals_x(jx) * vals_y(jy);
}
}
}
Expand Down
26 changes: 12 additions & 14 deletions Bembel/src/AnsatzSpace/SuperSpace.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,8 +15,8 @@ namespace Bembel {
* \ingroup AnsatzSpace
* \brief The superspace manages local polynomial bases on each element of the
* mesh and provides an interface to evaluate them.
*
*
*
*
*/
template <typename Derived>
struct SuperSpace {
Expand All @@ -26,7 +26,7 @@ struct SuperSpace {
//////////////////////////////////////////////////////////////////////////////
/**
* \brief Default constructor for the SuperSpace class.
*
*
* This constructor creates a SuperSpace object with default parameters.
*/
SuperSpace() {}
Expand Down Expand Up @@ -160,7 +160,8 @@ struct SuperSpace {
void map2surface(const ElementTreeNode& e, const Eigen::Vector2d& xi,
double w, SurfacePoint* surf_pt) const {
Eigen::Vector2d st = e.llc_ + e.get_h() * xi;
mesh_->get_geometry()[e.patch_].updateSurfacePoint(surf_pt, st, w, xi);
mesh_->get_geometry()[e.patch_].updateSurfacePoint(surf_pt, e.patch_, st, w,
xi);
return;
}
//////////////////////////////////////////////////////////////////////////////
Expand Down Expand Up @@ -196,23 +197,20 @@ struct SuperSpace {
void addScaledSurfaceCurlInteraction(
Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>* intval, Scalar w,
const SurfacePoint& p1, const SurfacePoint& p2) const {
// surface measures
double kappa1 = p1.segment<3>(6).cross(p1.segment<3>(9)).norm();
double kappa2 = p2.segment<3>(6).cross(p2.segment<3>(9)).norm();
// compute basis functions's surface curl. Each column of s_curl is a basis
// function's surface curl at point s.
Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> s_curl =
(1.0 / kappa1) *
(-p1.segment<3>(6) * basisDy(p1.segment<2>(0)).transpose() +
p1.segment<3>(9) * basisDx(p1.segment<2>(0)).transpose());
(1.0 / p1.get_surface_measure()) *
(-p1.get_f_dx() * basisDy(p1.get_xi()).transpose() +
p1.get_f_dy() * basisDx(p1.get_xi()).transpose());
Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> t_curl =
(1.0 / kappa2) *
(-p2.segment<3>(6) * basisDy(p2.segment<2>(0)).transpose() +
p2.segment<3>(9) * basisDx(p2.segment<2>(0)).transpose());
(1.0 / p2.get_surface_measure()) *
(-p2.get_f_dx() * basisDy(p2.get_xi()).transpose() +
p2.get_f_dy() * basisDx(p2.get_xi()).transpose());
// inner product of surface curls of any two basis functions
for (int j = 0; j < polynomial_degree_plus_one_squared; ++j)
for (int i = 0; i < polynomial_degree_plus_one_squared; ++i)
(*intval)(j * polynomial_degree_plus_one_squared + i) +=
(*intval)(j* polynomial_degree_plus_one_squared + i) +=
w * s_curl.col(i).dot(t_curl.col(j));
}

Expand Down
64 changes: 23 additions & 41 deletions Bembel/src/Geometry/Patch.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -140,7 +140,6 @@ class Patch {
const double tpbasisval = xbasis[i] * ybasis[j];
const int accs = 4 * (numy * (polynomial_degree_x_ * x_location + i) +
polynomial_degree_y_ * y_location + j);
#pragma omp simd
for (int k = 0; k < 4; k++) tmp[k] += data_[accs + k] * tpbasisval;
}
}
Expand Down Expand Up @@ -202,7 +201,6 @@ class Patch {

// Here I add up the values of the basis functions in the dc
// basis
#pragma omp simd
for (int k = 0; k < 4; k++) {
tmp[k] += data_[accs + k] * tpbasisval;
tmpDx[k] += data_[accs + k] * tpbasisvalDx;
Expand All @@ -216,13 +214,8 @@ class Patch {
delete[] xbasisD;
delete[] ybasisD;

Eigen::Matrix<double, 3, 2> out;

// Eigen::Vector3d out;

double bot = 1. / (tmp[3] * tmp[3]);

#pragma omp simd
Eigen::Matrix<double, 3, 2> out;
for (int k = 0; k < 3; k++) {
out(k, 0) = (tmpDx[k] * tmp[3] - tmp[k] * tmpDx[3]) * bot;
out(k, 1) = (tmpDy[k] * tmp[3] - tmp[k] * tmpDy[3]) * bot;
Expand Down Expand Up @@ -278,21 +271,10 @@ class Patch {
return evalNormal(Eigen::Vector2d(x, y));
}

//
/**
* \brief Updates the surface point and returns the physical point and the
* derivatives there.
*
* This is a combination of eval und evalJacobian, to avoid duplication of
* work.
*
* \param srf_pt Pointer to the SurfacePoint which gets updated.
* \param ref_pt Point in reference domain with respect to the patch.
* \param w quadrature weight.
* \param xi Point in reference domain with respect to the element.
*/
void updateSurfacePoint(SurfacePoint *srf_pt,
const Eigen::Vector2d &ref_pt, double w,
// This is a combination of eval und evalJacobian, to avoid duplication of
// work. See SurfacePoint.hpp
void updateSurfacePoint(SurfacePoint *srf_pt, const int patch,
const Eigen::Vector2d &ref_pt, const double w,
const Eigen::Vector2d &xi) const {
const int x_location =
Spl::FindLocationInKnotVector(ref_pt(0), unique_knots_x_);
Expand All @@ -304,11 +286,11 @@ class Patch {
const double scaledy = Spl::Rescale(ref_pt(1), unique_knots_y_[y_location],
unique_knots_y_[y_location + 1]);

// TODO(Felix) Do not use variable-length arrays in accordance to Google
// Style guide
double *buffer =
new double[2 * (polynomial_degree_x_ + polynomial_degree_y_) + 12];
for (int i = 0; i < 12; ++i) buffer[i] = 0;
// allocating memory is expensive, use the buffer in SurfacePoint
srf_pt->allocate_buffer(2 * (polynomial_degree_x_ + polynomial_degree_y_) +
12);
double *buffer = srf_pt->get_buffer();
std::memset(buffer, 0, 12 * sizeof(double));
Copy link

Copilot AI Jun 20, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Consider using the allocated buffer size rather than a fixed value (12 * sizeof(double)) in the memset call, so that all of the buffer is properly cleared when polynomial degrees change.

Suggested change
std::memset(buffer, 0, 12 * sizeof(double));
std::memset(buffer, 0, sizeof(double) * (2 * (polynomial_degree_x_ + polynomial_degree_y_) + 12));

Copilot uses AI. Check for mistakes.

double *tmp = buffer;
double *tmpDx = tmp + 4;
Expand Down Expand Up @@ -348,19 +330,19 @@ class Patch {
const double bot = 1. / tmp[3];
const double botsqr = bot * bot;

(*srf_pt)(0) = xi(0);
(*srf_pt)(1) = xi(1);
(*srf_pt)(2) = w;
(*srf_pt)(3) = tmp[0] * bot;
(*srf_pt)(4) = tmp[1] * bot;
(*srf_pt)(5) = tmp[2] * bot;
(*srf_pt)(6) = (tmpDx[0] * tmp[3] - tmp[0] * tmpDx[3]) * botsqr;
(*srf_pt)(7) = (tmpDx[1] * tmp[3] - tmp[1] * tmpDx[3]) * botsqr;
(*srf_pt)(8) = (tmpDx[2] * tmp[3] - tmp[2] * tmpDx[3]) * botsqr;
(*srf_pt)(9) = (tmpDy[0] * tmp[3] - tmp[0] * tmpDy[3]) * botsqr;
(*srf_pt)(10) = (tmpDy[1] * tmp[3] - tmp[1] * tmpDy[3]) * botsqr;
(*srf_pt)(11) = (tmpDy[2] * tmp[3] - tmp[2] * tmpDy[3]) * botsqr;
delete[] buffer;
Eigen::Vector3d tmpvec(tmp[0], tmp[1], tmp[2]);
Eigen::Vector3d tmpDxvec(tmpDx[0], tmpDx[1], tmpDx[2]);
Eigen::Vector3d tmpDyvec(tmpDy[0], tmpDy[1], tmpDy[2]);

srf_pt->set_patch(patch);
srf_pt->set_xi(xi);
srf_pt->set_xi_patch(ref_pt);
srf_pt->set_w(w);
srf_pt->set_f(bot * tmpvec);
srf_pt->set_jacobian(botsqr * (Eigen::Matrix<double, 3, 2>()
<< tmpDxvec * tmp[3] - tmpvec * tmpDx[3],
tmpDyvec * tmp[3] - tmpvec * tmpDy[3])
.finished());
return;
}

Expand Down
Loading