Commit b47c7aa1 authored by Jan Zapletal's avatar Jan Zapletal

ENH: updated readme, MAINT: new order of args in repr. formula

parent bcb3a0f5
......@@ -2,5 +2,8 @@ Examples/build/*
Examples/dist/*
build/*
dist/*
manual/formulas.*
manual/figs.svg
*.bak
.DS_Store
......@@ -114,13 +114,14 @@ void testHeatDtN(
double point [] = { 0.1, 0.2, 0.3 };
bem4i::evaluateLaplaceRepresentationFormula(
&value,
nNodes,
nodes,
nElems,
elems,
1,
point,
&value,
alpha,
dirichlet,
orderFar,
data,
......
......@@ -2,11 +2,13 @@
## About the program
The purpose of the library is to generate the Dirichlet-to-Neumann (Steklov-Poincaré) operator for the Laplace equation in 3D. It is primarily designed to be used in a boundary element tearing and interconnecting (BETI) solver The operator can be used to solve the stationary heat equation or an electrostatic problem and is primarily des. The operator is discretized by the Galerkin boundary element method.
The purpose of the library is to generate the Dirichlet-to-Neumann (Steklov-Poincaré) operator for the Laplace equation in 3D. The operator can be used to solve the stationary heat equation or an electrostatic problem and is primarily designed to be used in a boundary element tearing and interconnecting (BETI) solver. The operator is discretized by the Galerkin boundary element method.
## Mathematical background
We consider a Neumann boundary value problem for the Laplace equation in 3D,
### Local problem
In total (all-floating) BETI we consider a Neumann boundary value problem for the Laplace equation in 3D for every subdomain,
<img src="manual/bvp.png" alt="BVP"> with the compatibility condition <img src="manual/compat.png" alt="compatibility condition">.
......@@ -50,13 +52,27 @@ with the Neumann data discretized by piecewise constant functions and
<img src="manual/ai.png" alt="ai"/>, <img src="manual/a.png" alt="a"/>.
### Volume evaluation
By solving the global BETI problem we obtain the Dirichlet data. In order to evaluate the solution in every subdomain we solve the auxiliary Dirichlet problem
<img src="manual/bvpd.png" alt="BVP">
by the boundary integral equation
<img src="manual/dirbie.png" alt="BIE">
and use the representation formula
<img src="manual/repr.png" alt="representation formula">.
## Compilation
A Makefile to compile the project as a dynamic library `./dist/libheatdtn.so` is provided. The library relies on Eigen (tested with, e.g. version 3.3.6) and Intel MKL (or possibly another BLAS/LAPACK implementation) which are **not** distributed within this repository. The user is required to download both and set the `EIGEN_INC` variable in `./Makefile`, e.g.
```
EIGEN_INC=/home/zap150/eigen-eigen-b3f3d4950030
```
and run the MKL script `mklvars.sh` distributed with Intel MKL to se the `MKLROOT` variable. The variable `SIMD` defines the `simdlen` clause for OpenMP `simd` pragmas, i.e. should be set to `4` for AVX2 and `8` for AVX512. The variables `CXX`, `CXXFLAGS`, `LDLIBSOPTIONS` should be set by the user, see the examples below for the Intel and GNU toolchains. Please note that the Intel compiler is recommended to make use of vectorization.
and run the MKL script `mklvars.sh` distributed with Intel MKL to set the `MKLROOT` variable. The variable `SIMD` defines the `simdlen` argument for OpenMP `simd` pragmas, i.e. should be set to `4` for AVX2 and `8` for AVX512. The variables `CXX`, `CXXFLAGS`, `LDLIBSOPTIONS` should be set by the user, see the examples below for the Intel and GNU toolchains. Please note that the Intel compiler is recommended to make efficient use of vectorization.
### Intel
......@@ -118,24 +134,42 @@ void getLaplaceSteklovPoincare(
* `nonsymmetric` ... `false` for symmetric, `true` for non-symmetric version,
* `verbose` ... `false` for silent computation, `true` for progress info.
The single- and double- layer matrices are stored in `data` to be used for the evaluation of the representation formula.
### Representation formula
```
void evaluateLaplaceRepresentationFormula(
double * values,
int nNodes,
const double * nodes,
int nElems,
const int * elems,
int nPoints,
const double * points,
double * values,
double alpha,
double * dirichlet,
int qOrderDisj,
int orderFar,
bem4idata< int, double > *& data,
bool verbose = false
);
```
* `values` ... pointer to an **allocated** array of the size `nPoints`,
* `nNodes` ... number of mesh nodes,
* `nodes` ... pointer to an array of node coordinates `{ x1, y1, z1, x2, y2, z2, ... }`
* `nElems` ... number of triangular elements
* `elems` ... indices to the array of nodes defining triangular elements `{ i1, j1, k1, i2, j2, k2, ... }`, exterior normal direction to *i*-th element **must** agree with ((xj,yj,zj)-(xi,yi,zi))x((xk,yk,zk)-(xj,yj,zj)),
* `nPoints` ... number of evaluation points,
* `points` ... pointer to an array of node coordinates `{ x1, y1, z1, x2, y2, z2, ... }`,
* `alpha` ... heat conduction parameter,
* `dirichlet` ... pointer to an array containing the nodal values of the Dirichlet data,
* `orderFar` ... far-field quadrature order, choose from {0,...,7}, (`4` recommended),
* `data` ... pointer to heatDtN data,
* `verbose` ... `false` for silent computation, `true` for progress info.
The function can **only** be called after `getLaplaceSteklovPoincare`.
### Data destruction
```
......@@ -144,10 +178,19 @@ void deleteBem4iData(
);
```
The function has to be called to delete all heatDtN data correctly.
The function has to be called to delete all heatDtN data correctly. Note that `evaluateLaplaceRepresentationFormula` **cannot** be called after this call.
## heatDtN as an ESPRESO plugin
We refer to the documentation of the ESPRESO library at http://espreso.it4i.cz/, for further info please contact the developers of heatDtN or ESPRESO at jan.zapletal@vsb.cz or ondrej.meca@vsb.cz.
### Compilation
To include the heatDtN module include the parameter `HEATDTN::PATH = /path/to/heatdtn` in `build.config`.
## Compilation with ESPRESO
### Configuration
An example of a configuration file can be found [here](manual/heatdtn.ecf).
## Licence
......
......@@ -199,15 +199,16 @@ void getLaplaceSteklovPoincare(
template<class LO, class SC>
void evaluateLaplaceRepresentationFormula(
SC * values,
LO nNodes,
const SC * nodes,
LO nElems,
const LO * elems,
LO nPoints,
const SC * points,
SC * values,
SC alpha,
SC * dirichlet,
int qOrderDisj,
int orderFar,
bem4idata< LO, SC > *& data,
bool verbose
) {
......@@ -221,6 +222,7 @@ void evaluateLaplaceRepresentationFormula(
for ( LO i = 0; i < nPoints; ++i ) {
values[ i ] = 0.0;
}
return;
}
std::vector< SCVT > nodesv;
......@@ -242,7 +244,7 @@ void evaluateLaplaceRepresentationFormula(
V->CholeskiDecomposedSolve( neu );
RepresentationFormulaLaplace<LO, SC> formula( &bespace10, &dir, &neu,
qOrderDisj );
orderFar );
Vector< LO, SC > val( nPoints, values, false );
......@@ -275,15 +277,16 @@ template void bem4i::getLaplaceSteklovPoincare< int, double >(
);
template void bem4i::evaluateLaplaceRepresentationFormula< int, double >(
double * values,
int nNodes,
const double * nodes,
int nElems,
const int * elems,
int nPoints,
const double * points,
double * values,
double alpha,
double * dirichlet,
int qOrderDisj,
int orderFar,
bem4idata< int, double > *& data,
bool verbose
);
......@@ -314,15 +317,16 @@ template void bem4i::getLaplaceSteklovPoincare< long, double >(
);
template void bem4i::evaluateLaplaceRepresentationFormula< long, double >(
double * values,
long nNodes,
const double * nodes,
long nElems,
const long * elems,
long nPoints,
const double * points,
double * values,
double alpha,
double * dirichlet,
int qOrderDisj,
int orderFar,
bem4idata< int, double > *& data,
bool verbose
);
......
......@@ -9,7 +9,7 @@ template<class LO, class SC>
void deleteBem4iData(
bem4idata< LO, SC > *& data
);
template<class LO, class SC>
void getLaplaceSteklovPoincare(
SC * S,
......@@ -30,15 +30,16 @@ void getLaplaceSteklovPoincare(
template<class LO, class SC>
void evaluateLaplaceRepresentationFormula(
SC * values,
LO nNodes,
const SC * nodes,
LO nElems,
const LO * elems,
LO nPoints,
const SC * points,
SC * values,
SC alpha,
SC * dirichlet,
int qOrderDisj,
int orderFar,
bem4idata< LO, SC > *& data,
bool verbose = false
);
......@@ -46,4 +47,3 @@ void evaluateLaplaceRepresentationFormula(
}
#endif
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment