Commit 9ba74bb8 authored by Jan Zapletal's avatar Jan Zapletal

FIX: allow for no far-field, ENH: readme

parent 9183c7e8
......@@ -28,11 +28,11 @@ or in its non-symmetric form as
The individual matrices composing the operator are given as
<img src="manual/vh.png" alt="Vh"/>,
<img src="manual/vh.png" alt="Vh"/>,
<img src="manual/kh.png" alt="Kh"/>,
<img src="manual/kh.png" alt="Kh"/>,
<img src="manual/dh.png" alt="Dh"/>,
<img src="manual/dh.png" alt="Dh"/>,
<img src="manual/mh.png" alt="Mh"/>
......@@ -44,14 +44,12 @@ can be obtained from the discretized boundary integral equation
<img src="manual/bem.png" alt="BEM"/>
with
with the Neumann data discretized by piecewise constant functions and
<img src="manual/beta.png" alt="beta"/>,
<img src="manual/beta.png" alt="beta"/>,
<img src="manual/ai.png" alt="ai"/>, <img src="manual/a.png" alt="a"/>.
## 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.
......@@ -82,6 +80,8 @@ A simple program linking against `libheatdtn.so` is provided in the directory `.
## heatDtN API
The library provides three functions described below.
### Generation of the Dirichlet-to-Neumann map
```
......@@ -103,6 +103,21 @@ void getLaplaceSteklovPoincare(
);
```
* `S` ... pointer to an **allocated** array of the size `nNodes * nNodes`,
* `a` ... pointer to an **allocated** array of the size `nNodes`,
* `beta` ... pointer to an **allocated** array of the size `1`,
* `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))
* `alpha` ... heat conduction parameter,
* `qType` ... `0` for semi-analytic quadrature, `1` for regularized numerical scheme,
* `orderNear` ... near-field quadrature order, choose from {1,...,7} for `qType = 0` (`5` recommended), {3,...,10} for `qType = 1` (`4` recommended),
* `orderFar` ... far-field quadrature order, choose from {0,...,7}, (`4` recommended, `0` uses near-field quadrature for far-field as well),
* `data` ... pointer to heatDtN data,
* `nonsymmetric` ... `false` for symmetric, `true` for non-symmetric version,
* `verbose` ... `false` for silent computation, `true` for progress info.
### Representation formula
```
......@@ -121,7 +136,13 @@ void evaluateLaplaceRepresentationFormula(
);
```
### Data destruction
```
void deleteBem4iData(
bem4idata< LO, SC > *& data
);
```
## Compilation with ESPRESO
......
......@@ -59,7 +59,7 @@ void getLaplaceSteklovPoincare(
const SC * nodes,
LO nElems,
const LO * elems,
SC alpha,
SC alpha,
int qType,
int orderNear,
int orderFar,
......@@ -69,13 +69,13 @@ void getLaplaceSteklovPoincare(
) {
typedef typename GetType<LO, SC>::SCVT SCVT;
if ( qType < 0 || qType > 1 ) qType = 1;
if ( qType == 1 && orderNear < 1 ) orderNear = 1;
if ( qType == 1 && orderNear > 7 ) orderNear = 7;
if ( qType == 0 && orderNear < 3 ) orderNear = 3;
if ( qType == 0 && orderNear > 10 ) orderNear = 10;
if ( orderFar < 1 ) orderFar = 1;
if ( orderFar < 0 ) orderFar = 0;
if ( orderFar > 7 ) orderFar = 7;
data = new bem4idata< LO, SC >( );
......@@ -141,15 +141,15 @@ void getLaplaceSteklovPoincare(
Vector< LO, SC > onese( nElems );
onese.setAll( 1.0 );
M.apply( onese, avect, true, 1.0, 0.0 );
if( beta != nullptr ){
Vector< LO, SC > onesn( nNodes );
onesn.setAll( 1.0 );
*beta = avect.dot( onesn );
*beta = 0.25 / *beta;
*beta = 0.25 / *beta;
}
}
Eigen::SparseMatrix<SC, Eigen::ColMajor, LO> * Me = M.getEigenSparseMatrix( );
// K <- 1/2M + K
......@@ -238,7 +238,7 @@ void evaluateLaplaceRepresentationFormula(
// rhs = (1/2 I + K)*dir
K->apply( dir, neu, false, 1.0, 0.0 );
// local solve (V decomposed in assembly)
// local solve (V decomposed in assembly)
V->CholeskiDecomposedSolve( neu );
RepresentationFormulaLaplace<LO, SC> formula( &bespace10, &dir, &neu,
......@@ -265,7 +265,7 @@ template void bem4i::getLaplaceSteklovPoincare< int, double >(
const double * nodes,
int nElems,
const int * elems,
double alpha,
double alpha,
int qType,
int orderNear,
int orderFar,
......@@ -328,4 +328,3 @@ template void bem4i::evaluateLaplaceRepresentationFormula< long, double >(
);
#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