Skip to content
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.

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

Remove auto variables in tangential trace #70

Merged
merged 7 commits into from
Feb 13, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
1 change: 1 addition & 0 deletions Bembel/LinearForm
Original file line number Diff line number Diff line change
Expand Up @@ -26,5 +26,6 @@
#include "src/LinearForm/DirichletTrace.hpp"
#include "src/LinearForm/DiscreteLinearForm.hpp"
#include "src/LinearForm/RotatedTangentialTrace.hpp"
#include "src/LinearForm/TangentialTrace.hpp"

#endif // BEMBEL_LINEARFORM_MODULE_
88 changes: 88 additions & 0 deletions Bembel/src/LinearForm/TangentialTrace.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,88 @@
// This file is part of Bembel, the higher order C++ boundary element library.
//
// Copyright (C) 2022 see <http://www.bembel.eu>
//
// It was written as part of a cooperation of J. Doelz, H. Harbrecht, S. Kurz,
// M. Multerer, S. Schoeps, and F. Wolf at Technische Universitaet Darmstadt,
// Universitaet Basel, and Universita della Svizzera italiana, Lugano. This
// source code is subject to the GNU General Public License version 3 and
// provided WITHOUT ANY WARRANTY, see <http://www.bembel.eu> for further
// information.
#ifndef BEMBEL_SRC_LINEARFORM_TANGENTIALTRACE_HPP_
#define BEMBEL_SRC_LINEARFORM_TANGENTIALTRACE_HPP_

namespace Bembel {

template <typename Scalar>
class TangentialTrace;

template <typename ScalarT>
struct LinearFormTraits<TangentialTrace<ScalarT>> {
typedef ScalarT Scalar;
};

/**
* \ingroup LinearForm
* \brief This class provides a specialization of the linear form required
*for the solution of the electric field integral equation.
**/
template <typename Scalar>
class TangentialTrace : public LinearFormBase<TangentialTrace<Scalar>, Scalar> {
public:
TangentialTrace() {}
void set_function(
const std::function<Eigen::Matrix<Scalar, 3, 1>(Eigen::Vector3d)>
&function) {
function_ = function;
}
template <class T>
void evaluateIntegrand_impl(
const T &super_space, const SurfacePoint &p,
Eigen::Matrix<Scalar, Eigen::Dynamic, 2> *intval) const {
int polynomial_degree = super_space.get_polynomial_degree();
int polynomial_degree_plus_one_squared =
(polynomial_degree + 1) * (polynomial_degree + 1);

// get evaluation points on unit square
Eigen::Vector2d s = p.segment<2>(0);

// get quadrature weights
double ws = p(2);

// get points on geometry and tangential derivatives
Eigen::Vector3d x_f = p.segment<3>(3);
Eigen::Vector3d x_f_dx = p.segment<3>(6);
Eigen::Vector3d x_f_dy = p.segment<3>(9);

// compute surface measures from tangential derivatives
Eigen::Vector3d x_n = x_f_dx.cross(x_f_dy).normalized();

// tangential component + quadrature weights
Eigen::Matrix<Scalar, 3, 1> fun_x_f = function_(x_f);
Eigen::Matrix<Scalar, 3, 1> tangential_component = fun_x_f.cross(x_n) * ws;

// extract tangential component
Scalar component_x = x_f_dx.dot(tangential_component);
Scalar component_y = x_f_dy.dot(tangential_component);

// evaluate shape functions
Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> phiPhiVec =
super_space.basis(s);

// multiply basis functions with integrand
Eigen::Matrix<Scalar, Eigen::Dynamic, 2> phiPhiMat(
polynomial_degree_plus_one_squared, 2);
phiPhiMat.col(0) = component_x * phiPhiVec;
phiPhiMat.col(1) = component_y * phiPhiVec;

// compute integrals
(*intval) += phiPhiMat;
return;
}

private:
std::function<Eigen::Matrix<Scalar, 3, 1>(Eigen::Vector3d)> function_;
};
} // namespace Bembel

#endif // BEMBEL_SRC_LINEARFORM_TANGENTIALTRACE_HPP_
1 change: 1 addition & 0 deletions tests/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@ set(UNITTESTS
test_SurfacePointUpdate
test_Projector
test_Glue
test_TangentialTrace
test_FMMTransferMatrices
test_FMMForwardTransformation
test_FMMBackwardTransformation
Expand Down
51 changes: 51 additions & 0 deletions tests/test_TangentialTrace.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,51 @@
// This file is part of Bembel, the higher order C++ boundary element library.
//
// Copyright (C) 2024 see <http://www.bembel.eu>
//
// It was written as part of a cooperation of J. Doelz, H. Harbrecht, S. Kurz,
// M. Multerer, S. Schoeps, and F. Wolf at Technische Universitaet Darmstadt,
// Universitaet Basel, and Universita della Svizzera italiana, Lugano. This
// source code is subject to the GNU General Public License version 3 and
// provided WITHOUT ANY WARRANTY, see <http://www.bembel.eu> for further
// information.

#include <Bembel/AnsatzSpace>
#include <Bembel/Geometry>
#include <Bembel/LinearForm>

#include "tests/TestGeometries.hpp"

class TestOperatorDivC;
template <>
struct Bembel::LinearOperatorTraits<TestOperatorDivC> {
typedef Eigen::VectorXd EigenType;
typedef Eigen::VectorXd::Scalar Scalar;
enum { OperatorOrder = 0, Form = DifferentialForm::DivConforming };
};

int main() {
Test::TestGeometryWriter::writeScreen();
Bembel::Geometry geometry("test_Screen.dat");

const int refinement_level = 0;
const int polynomial_degree = 1;
Bembel::AnsatzSpace<TestOperatorDivC> ansatz_space(geometry, refinement_level,
polynomial_degree);

std::function<Eigen::Vector3d(Eigen::Vector3d)> fun = [](Eigen::Vector3d in) {
return Eigen::Vector3d(in(0), in(1), in(2));
};

Bembel::DiscreteLinearForm<Bembel::TangentialTrace<double>, TestOperatorDivC>
disc_lf(ansatz_space);
disc_lf.get_linear_form().set_function(fun);
disc_lf.compute();

Eigen::VectorXd ref_sol(4);
ref_sol << 0.25, 0.25, -0.25, -0.25;

BEMBEL_TEST_IF((disc_lf.get_discrete_linear_form() - ref_sol).norm() <
Test::Constants::coefficient_accuracy);

return 0;
}
Loading