Skip to content

Commit

Permalink
Add Augmented Electric Field Integral Equation Feature (#101)
Browse files Browse the repository at this point in the history
* Add Screen Feature

* Add AugmentedEFIE Feature

* Add AugmentedEFIE to Documentation

* Add Readme to Unsupported Directory

* Install Specific cpplint Version Without Namespace Bug
  • Loading branch information
mx-nlte authored Dec 18, 2024
1 parent 82e998a commit 1e21364
Show file tree
Hide file tree
Showing 26 changed files with 4,219 additions and 15 deletions.
2 changes: 1 addition & 1 deletion .github/workflows/code-style.yml
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ jobs:
- name: Install dependencies
run: |
python -m pip install --upgrade pip
pip install cpplint
pip install --force-reinstall -v "cpplint==1.6.1"
- name: "Check Code Format"
run: |
./scripts/cpplint.sh
33 changes: 25 additions & 8 deletions Bembel/src/AnsatzSpace/Glue.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -165,13 +165,15 @@ class Glue {

// Here, we fill the two vectors above
int number_of_slaves = 0;
int number_of_boundary_dofs = 0;
for (auto dofset : dof_id) {
// This sorting is required, since by construction only dofs[0]<dofs[1] is
// given, but in the case of corners and H1 discretizations, it might
// happen that dofs[0]>dofs[2]. Moreover, we count how many dofs we "glue
// away".
// away" or skipped on the boundary.
std::sort(dofset.dofs.begin(), dofset.dofs.end());
dof_is_master[dofset.dofs[0]] = true;
if (dofset.dofs.size() == 1) ++number_of_boundary_dofs;
for (int i = 1; i < dofset.dofs.size(); ++i) {
dof_is_slave[dofset.dofs[i]] = true;
++number_of_slaves;
Expand All @@ -185,17 +187,24 @@ class Glue {
return a.dofs[0] < b.dofs[0];
});

const int post_dofs = pre_dofs - number_of_slaves;

const int post_dofs = pre_dofs - number_of_slaves - number_of_boundary_dofs;
const int glued_dofs = pre_dofs - number_of_slaves;
assert(post_dofs != 0 && "All degrees of freedom are on the boundary!");
// This block is the heart of the algorithms. Skip keeps track of how many
// slaves have been skipped already, and master_index keeps track on which
// master is about to be glued next.
// master is about to be glued next. post_index keeps track on how many
// dofs are on the boundary of the patches and therefore skipped.
int skip = 0;
int master_index = 0;
for (int i = 0; i < post_dofs; ++i) {
const int post_index = i;
int post_index = 0;
for (int i = 0; i < glued_dofs; ++i) {
while (dof_is_slave[i + skip] && ((i + skip) < pre_dofs)) ++skip;
const int pre_index = i + skip;
// skip dofs on patch boundary without a slave patch
if (dof_is_master[pre_index] && dof_id[master_index].dofs.size() == 1) {
++master_index;
continue;
}
trips.push_back(Eigen::Triplet<double>(pre_index, post_index, 1));
// The dof cannot be declared slave and master at the same time
assert(!(dof_is_master[pre_index] && dof_is_slave[pre_index]));
Expand All @@ -212,12 +221,13 @@ class Glue {
}
++master_index;
}
++post_index;
}

Eigen::SparseMatrix<double> glue_matrix(pre_dofs, post_dofs);
glue_matrix.setFromTriplets(trips.begin(), trips.end());

assert(trips.size() == pre_dofs);
assert(trips.size() == (pre_dofs - number_of_boundary_dofs));
return glue_matrix;
}
};
Expand Down Expand Up @@ -557,7 +567,7 @@ struct glue_identificationmaker_<Derived, DifferentialForm::DivConforming> {
const bool needReversion = reverseParametrized(edge);
const int coef = glueCoefficientDivergenceConforming(edge);

// Check if the edge is hanging. For screens one would need an "else".
// Check if the edge is hanging.
if (edge[1] > -1 && edge[0] > -1) {
for (int i = 0; i < size_of_edge_dofs; ++i) {
GlueRoutines::dofIdentification d;
Expand All @@ -569,6 +579,13 @@ struct glue_identificationmaker_<Derived, DifferentialForm::DivConforming> {
d.coef = coef;
out.push_back(d);
}
} else {
for (int i = 0; i < size_of_edge_dofs; ++i) {
// Add only the master dof which is skipped in subsequent routines
GlueRoutines::dofIdentification d;
d.dofs.push_back(dofs_e1[i]);
out.push_back(d);
}
}
}
out.shrink_to_fit();
Expand Down
5 changes: 5 additions & 0 deletions Bembel/src/util/Constants.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,11 @@ constexpr size_t Bembel_alloc_size = 100;
// solution of the linear system
constexpr double projector_tolerance = 1e-4;
////////////////////////////////////////////////////////////////////////////////
/// physical constants
////////////////////////////////////////////////////////////////////////////////
constexpr double eps0 = 8.854187812813 * 1e-12;
constexpr double mu0 = 4 * M_PI * 1e-7;
////////////////////////////////////////////////////////////////////////////////
/// methods
////////////////////////////////////////////////////////////////////////////////

Expand Down
2 changes: 2 additions & 0 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -73,6 +73,7 @@ set (BEMBEL_PATH "${PROJECT_SOURCE_DIR}")

set(CMAKE_CXX_STANDARD 11)
include_directories(${BEMBEL_PATH})
include_directories(${BEMBEL_PATH}/unsupported)

if(NOT CMAKE_BUILD_TYPE)
set(CMAKE_BUILD_TYPE Release CACHE STRING
Expand Down Expand Up @@ -108,4 +109,5 @@ set(CMAKE_CXX_FLAGS_COVERAGE "-g -fPIC -fprofile-arcs -ftest-coverage" CACHE STR

add_subdirectory(tests)
add_subdirectory(examples)
add_subdirectory(unsupported/examples)
enable_testing()
3 changes: 2 additions & 1 deletion Doxyfile
Original file line number Diff line number Diff line change
Expand Up @@ -833,7 +833,8 @@ WARN_LOGFILE =
INPUT = Bembel \
README.md \
doc \
examples
examples \
unsupported

# This tag can be used to specify the character encoding of the source files
# that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses
Expand Down
2 changes: 1 addition & 1 deletion doc/layout.xml
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@
<tab type="user" url="citelist.html" title="Publications Featuring Bembel"/>
<tab type="user" url="@ref Funding" title="Funding"/>
</tab>
<tab type="modules" visible="yes" title="Source Code" intro=""/>
<tab type="topics" visible="yes" title="Source Code" intro=""/>
<tab type="classlist" visible="yes" title="" intro=""/>
<tab type="namespacelist" visible="yes" title=""/>
</navindex>
Expand Down
31 changes: 27 additions & 4 deletions tests/test_TangentialTrace.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -22,18 +22,41 @@ struct Bembel::LinearOperatorTraits<TestOperatorDivC> {
typedef Eigen::VectorXd::Scalar Scalar;
enum { OperatorOrder = 0, Form = DifferentialForm::DivConforming };
};

/*
* This test uses a five patch geometry. After discretizing with level zero
* refinement four degrees of freedom remain which are all located on the edge
* of patch 0.
* The linear form integrates the function (x,y,z) with the tangential trace.
* For a simplified test the function is 0 outside of patch 0. Thereby only
* contributions of basis functions with support on patch 0 need to be taken
* into account.
*
* ----
* |4 |
* ----------
* |2 |0 |1 |
* ----------
* |3 |
* ----
*
*/
int main() {
Test::TestGeometryWriter::writeScreen();
Bembel::Geometry geometry("test_Screen.dat");
Test::TestGeometryWriter::writeEdgeCase1();
Bembel::Geometry geometry("test_EdgeCase1.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));
Eigen::Vector3d retval;
if (in(0) < 0 || in(0) > 1 || in(1) < 0 || in(1) > 1) {
retval << 0, 0, 0;
} else {
retval << in(0), in(1), in(2);
}
return retval;
};

Bembel::DiscreteLinearForm<Bembel::TangentialTrace<double>, TestOperatorDivC>
Expand Down
38 changes: 38 additions & 0 deletions unsupported/Bembel/AugmentedEFIE
Original file line number Diff line number Diff line change
@@ -0,0 +1,38 @@
// 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.

#ifndef BEMBEL_AugmentedEFIE_MODULE_
#define BEMBEL_AugmentedEFIE_MODULE_

/**
* \ingroup Modules
* \defgroup AugmentedEFIE AugmentedEFIE
* \brief The AugmentedEFIE module contains routines for the spline-based
* Augmented Electric Field Integral Equation.
**/

#include <Bembel/LinearOperator>
#include <Bembel/Identity>
#include <Bembel/LinearForm>
#include <Bembel/Helmholtz>

#include "src/AugmentedEFIE/IncidenceMatrix.hpp"
#include "src/AugmentedEFIE/InductanceMatrix.hpp"
#include "src/AugmentedEFIE/EFIE.hpp"
#include "src/AugmentedEFIE/VectorPotential.hpp"
#include "src/AugmentedEFIE/GradScalarPotential.hpp"
#include "src/AugmentedEFIE/MixedEFIE.hpp"
#include "src/AugmentedEFIE/VoltageSource.hpp"
#include "src/AugmentedEFIE/AugmentedEFIEAssembler.hpp"
#include "src/AugmentedEFIE/AugmentedEFIE.hpp"
#include "src/AugmentedEFIE/AugmentedEFIEExcitation.hpp"

#endif
167 changes: 167 additions & 0 deletions unsupported/Bembel/src/AugmentedEFIE/AugmentedEFIE.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,167 @@
// 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.
#ifndef BEMBEL_SRC_AUGMENTEDEFIE_AUGMENTEDEFIE_HPP_
#define BEMBEL_SRC_AUGMENTEDEFIE_AUGMENTEDEFIE_HPP_

namespace Bembel {
/**
* \ingroup AugmentedEFIE
* \brief This class handles the assembly and storage of the system matrix.
*/
template <typename MatrixFormat, typename Derived>
class AugmentedEFIE {
public:
//////////////////////////////////////////////////////////////////////////////
// constructors
//////////////////////////////////////////////////////////////////////////////
AugmentedEFIE() {}
AugmentedEFIE(const AnsatzSpace<Derived> &ansatz_space_vector,
const Geometry &geometry) {
init_AugmentedEFIE(ansatz_space_vector, geometry);
}
//////////////////////////////////////////////////////////////////////////////
// init_AugmentedEFIE
//////////////////////////////////////////////////////////////////////////////
void init_AugmentedEFIE(const AnsatzSpace<Derived> &ansatz_space_vector,
const Geometry &geometry) {
ansatz_space_vector_ = ansatz_space_vector;

// Build ansatz spaces
ansatz_space_mass_ = AnsatzSpace<MassMatrixScalarDisc>(
geometry, ansatz_space_vector.get_refinement_level(),
ansatz_space_vector.get_polynomial_degree() - 1);

ansatz_space_scalar_ = AnsatzSpace<HelmholtzSingleLayerOperator>(
geometry, ansatz_space_vector.get_refinement_level(),
ansatz_space_vector.get_polynomial_degree() - 1);

dofs_scalar_ = ansatz_space_scalar_.get_number_of_dofs();
dofs_vector_ = ansatz_space_vector_.get_number_of_dofs();
return;
}
//////////////////////////////////////////////////////////////////////////////
// compute
//////////////////////////////////////////////////////////////////////////////
/**
* \brief Assembles the system matrix without voltage source excitation.
*/
inline void compute() {
system_matrix_ = Eigen::MatrixXcd(dofs_scalar_ + dofs_vector_,
dofs_scalar_ + dofs_vector_);

AugmentedEFIEAssembler<MatrixFormat>::compute(
&system_matrix_, ansatz_space_mass_, ansatz_space_scalar_,
ansatz_space_vector_, wavenumber_);
return;
}
/**
* \brief Assembles the system matrix with voltage source excitation.
*/
inline void compute(VoltageSource source) {
const int elements_per_edge =
(1 << ansatz_space_vector_.get_refinement_level());
// This number contains all degrees of freedom on the edge for the
// excitation and three extra ones dedicated to the potential of the voltage
// source.
const int dofs_excitation =
3 +
(source.positive_.size() + source.negative_.size()) * elements_per_edge;
system_matrix_ =
Eigen::MatrixXcd::Zero(dofs_scalar_ + dofs_vector_ + dofs_excitation,
dofs_scalar_ + dofs_vector_ + dofs_excitation);
excitation_ = Eigen::VectorXcd::Zero(dofs_scalar_);

AugmentedEFIEAssembler<MatrixFormat>::compute_matrix_and_excitation(
&system_matrix_, &excitation_, ansatz_space_mass_, ansatz_space_scalar_,
ansatz_space_vector_, wavenumber_, source);
return;
}
void stabilize() {
std::function<std::complex<double>(Eigen::Vector3d)> one =
[](Eigen::Vector3d in) { return std::complex<double>(1., 0.); };
DiscreteLinearForm<DirichletTrace<std::complex<double>>,
HelmholtzSingleLayerOperator>
const_lf(ansatz_space_scalar_);
const_lf.get_linear_form().set_function(one);
const_lf.compute();

system_matrix_.bottomRows(dofs_scalar_) *= omega_ *
std::complex<double>(0., 1.) *
Constants::mu0 * Constants::eps0;
system_matrix_.leftCols(dofs_vector_) *=
1. / (omega_ * std::complex<double>(0., 1.) * Constants::mu0);

std::complex<double> eigenvalues_average =
system_matrix_.trace() / ((double)system_matrix_.cols());

Eigen::VectorXcd a = Eigen::VectorXcd::Zero(system_matrix_.cols());
a.bottomRows(dofs_scalar_) = const_lf.get_discrete_linear_form();

system_matrix_ = system_matrix_ - eigenvalues_average * a * a.transpose();
}
//////////////////////////////////////////////////////////////////////////////
// setters
//////////////////////////////////////////////////////////////////////////////
/**
* \brief Set wave number
*/
void set_wavenumber(std::complex<double> wavenumber) {
wavenumber_ = wavenumber;
}
void set_omega(double omega) { omega_ = omega; }
//////////////////////////////////////////////////////////////////////////////
// getters
//////////////////////////////////////////////////////////////////////////////
const MatrixFormat &get_system_matrix() const { return system_matrix_; }
/**
* \brief Return wave number
*/
const std::complex<double> get_wavenumber() { return wavenumber_; }
/**
* \brief Return number of degrees of freedom for the current.
*/
const int get_dofs_vector() { return dofs_vector_; }
/**
* \brief Return number of degrees of freedom for the potential
*/
const int get_dofs_scalar() { return dofs_scalar_; }
/**
* \brief Debug output for the excitation.
*/
const Eigen::VectorXcd get_excitation() { return excitation_; }
/**
* \brief Debug output for the excitation.
*/
const AnsatzSpace<HelmholtzSingleLayerOperator> get_ansatz_space_scalar() {
return ansatz_space_scalar_;
}
const AnsatzSpace<MassMatrixScalarDisc> get_ansatz_space_mass() {
return ansatz_space_mass_;
}
//////////////////////////////////////////////////////////////////////////////
// private member variables
//////////////////////////////////////////////////////////////////////////////
private:
MatrixFormat system_matrix_;
Eigen::VectorXcd excitation_;

AnsatzSpace<MassMatrixScalarDisc> ansatz_space_mass_;
AnsatzSpace<HelmholtzSingleLayerOperator> ansatz_space_scalar_;
AnsatzSpace<InductanceMatrix> ansatz_space_vector_;

int dofs_scalar_;
int dofs_vector_;
double omega_;
std::complex<double> wavenumber_;
};

} // namespace Bembel
#endif // BEMBEL_SRC_AUGMENTEDEFIE_AUGMENTEDEFIE_HPP_
Loading

0 comments on commit 1e21364

Please sign in to comment.