Commit 7985d943 authored by Jan Zapletal's avatar Jan Zapletal

initial working commit

parent 5b37079e
Examples/testHeatDtN
build/*
dist/*
.dep.inc
.DS_Store
*.o
/*!
* @file BEBilinearForm.cpp
* @author Michal Merta
* @date July 12, 2013
*
*/
#ifdef BEBILINEARFORM_H
namespace bem4i {
template<class LO, class SC>
BEBilinearForm<LO, SC>::BEBilinearForm( ) {
}
template<class LO, class SC>
BEBilinearForm<LO, SC>::BEBilinearForm(
const BEBilinearForm& orig
) {
}
template<class LO, class SC>
BEBilinearForm<LO, SC>::~BEBilinearForm( ) {
}
}
#endif
/*!
* @file BEBilinearForm.h
* @author Michal Merta
* @date July 12, 2013
* @brief Header file for pure virtual class BEBilinearForm
*
*/
#ifndef BEBILINEARFORM_H
#define BEBILINEARFORM_H
#include <algorithm>
#include <omp.h>
#include <new>
#include "FullMatrix.h"
#include "Settings.h"
#include "Macros.h"
#include "Quadratures.h"
#include "ProgressMonitor.h"
namespace bem4i {
//! type of quadrature
enum quadratureType {
SauterSchwab,
Steinbach
};
/*!
* Abstract class representing a bilinear form over some BEM space
*
*/
template<class LO, class SC>
class BEBilinearForm {
typedef typename GetType<LO, SC>::SCVT SCVT;
public:
//! default constructor
BEBilinearForm( );
//! copy constructor
BEBilinearForm(
const BEBilinearForm& orig
);
//! destructor
virtual ~BEBilinearForm( );
//! method assembles a full matrix of a given bilinear form
virtual void assemble(
FullMatrix<LO, SC>& matrix
) const = 0;
protected:
//! boundary element space for full BEM
BESpace<LO, SC>* space;
//! quadrature rule
int* quadratureOrder;
//! quadrature type
quadratureType quadrature;
//! disjoint elements Gaussian quadrature order
int* quadratureOrderDisjointElems;
};
}
// include .cpp file to overcome linking problems due to templates
#include "BEBilinearForm.cpp"
#endif /* BEBILINEARFORM_H */
/*!
* @file BEBilinearFormLaplace1Layer.cpp
* @author Michal Merta
* @date August 8, 2013
*
*/
#ifdef BEBILINEARFORMLAPLACE1LAYER_H
namespace bem4i {
template<class LO, class SC>
BEBilinearFormLaplace1Layer<LO, SC>::BEBilinearFormLaplace1Layer( ) {
}
template<class LO, class SC>
BEBilinearFormLaplace1Layer<LO, SC>::BEBilinearFormLaplace1Layer(
const BEBilinearFormLaplace1Layer& orig
) {
this->space = orig.space;
this->quadratureOrder = orig.quadratureOrder;
this->quadratureOrderDisjointElems = orig.quadratureOrderDisjointElems;
}
template<class LO, class SC>
BEBilinearFormLaplace1Layer<LO, SC>::BEBilinearFormLaplace1Layer(
BESpace<LO, SC>* space,
int* quadratureOrder,
quadratureType quadrature,
int* quadratureOrderDisjointElems,
bool employSymmetricity
) {
this->space = space;
this->quadrature = quadrature;
if ( quadratureOrder ) {
this->quadratureOrder = quadratureOrder;
} else {
switch ( quadrature ) {
case SauterSchwab:
this->quadratureOrder = defaultQuadraturesSauterSchwab;
break;
case Steinbach:
this->quadratureOrder = defaultQuadraturesSteinbach;
}
}
this->quadratureOrderDisjointElems = quadratureOrderDisjointElems;
}
template<class LO, class SC>
BEBilinearFormLaplace1Layer<LO, SC>::~BEBilinearFormLaplace1Layer( ) {
}
template<class LO, class SC>
void BEBilinearFormLaplace1Layer<LO, SC>::assemble(
FullMatrix<LO, SC> & matrix
) const {
if ( this->space->getTestFunctionType( ) == p0 &&
this->space->getAnsatzFunctionType( ) == p0 ) {
this->assembleP0P0( matrix );
return;
} else if ( this->space->getTestFunctionType( ) == p1 &&
this->space->getAnsatzFunctionType( ) == p1 ) {
this->assembleP1P1( matrix );
return;
} else if ( this->space->getTestFunctionType( ) == p1dis &&
this->space->getAnsatzFunctionType( ) == p1dis ) {
this->assembleP1DisP1Dis( matrix );
return;
}
}
template<class LO, class SC>
void BEBilinearFormLaplace1Layer<LO, SC>::assembleP0P0(
FullMatrix<LO, SC> & matrix
) const {
LO nElems = this->space->getMesh( )->getNElements( );
matrix.resize( nElems, nElems, false );
SC * matrixData = matrix.getData( );
#pragma omp parallel
{
BEIntegratorLaplace<LO, SC> integrator( this->space, this->quadratureOrder,
this->quadrature, this->quadratureOrderDisjointElems );
FullMatrix< LO, SC > elemMatrix( 1, 1 );
#pragma omp for schedule( dynamic, 8 )
for ( LO j = 0; j < nElems; ++j ) {
for ( LO i = 0; i < nElems; ++i ) {
integrator.getElemMatrix1Layer( i, j, elemMatrix );
matrixData[ j * nElems + i ] = elemMatrix.get( 0, 0 );
}
}
} // end omp parallel
}
template<class LO, class SC>
void BEBilinearFormLaplace1Layer<LO, SC>::assembleP1P1(
FullMatrix<LO, SC> & matrix
) const {
LO nNodes = this->space->getMesh( )->getNNodes( );
LO nElems = this->space->getMesh( )->getNElements( );
matrix.resize( nNodes, nNodes, false );
SC * matrixData = matrix.getData( );
#pragma omp parallel for schedule( static, 32 )
for ( LO j = 0; j < nNodes; ++j ) {
for ( LO i = 0; i < nNodes; ++i ) {
matrixData[ j * nNodes + i ] = 0.0;
}
}
#pragma omp parallel
{
BEIntegratorLaplace<LO, SC> integrator( this->space,
this->quadratureOrder, this->quadrature,
this->quadratureOrderDisjointElems );
FullMatrix< LO, SC > elemMatrix( 3, 3 );
LO iElem[ 3 ], oElem[ 3 ];
SC * elemMatrixData = elemMatrix.getData( );
#pragma omp for schedule( dynamic, 8 )
for ( LO i = 0; i < nElems; ++i ) { // i outer row
for ( LO j = 0; j < nElems; ++j ) { // j inner col
integrator.getElemMatrix1Layer( i, j, elemMatrix );
this->space->getMesh( )->getElement( i, oElem );
this->space->getMesh( )->getElement( j, iElem );
for ( int oRot = 0; oRot < 3; ++oRot ) {
for ( int iRot = 0; iRot < 3; ++iRot ) {
#pragma omp atomic update
matrixData[ iElem[ iRot ] * nNodes + oElem[ oRot ] ] +=
elemMatrixData[ iRot * 3 + oRot ];
}
}
}
}
} // end omp parallel
}
template<class LO, class SC>
void BEBilinearFormLaplace1Layer<LO, SC>::assembleP1DisP1Dis(
FullMatrix<LO, SC> & matrix
) const {
LO nElems = this->space->getMesh( )->getNElements( );
matrix.resize( 3 * nElems, 3 * nElems, false );
SC * matrixData = matrix.getData( );
#pragma omp parallel
{
BEIntegratorLaplace<LO, SC> integrator( this->space,
this->quadratureOrder, this->quadrature,
this->quadratureOrderDisjointElems );
FullMatrix< LO, SC > elemMatrix( 3, 3 );
SC * elemMatrixData = elemMatrix.getData( );
#pragma omp for //schedule( dynamic, 8 )
for ( LO j = 0; j < nElems; ++j ) { // j inner col
for ( LO i = 0; i < nElems; ++i ) { // i outer row
integrator.getElemMatrix1Layer( i, j, elemMatrix );
for ( int oRot = 0; oRot < 3; ++oRot ) {
for ( int iRot = 0; iRot < 3; ++iRot ) {
matrixData[ ( 3 * j + iRot ) * 3 * nElems + 3 * i + oRot ] =
elemMatrixData[ iRot * 3 + oRot ];
}
}
}
}
} // end omp parallel
}
}
#endif
/*!
* @file BEBilinearFormLaplace1Layer.h
* @author Michal Merta
* @date July 12, 2013
* @brief Header file for class BEBilinearFormLaplace1Layer
*
*/
#ifndef BEBILINEARFORMLAPLACE1LAYER_H
#define BEBILINEARFORMLAPLACE1LAYER_H
#include "BEBilinearForm.h"
namespace bem4i {
/*!
* Class representing the bilinear form for the Laplace single layer operator
*
* Provides methods for system matrix assembly.
*
*/
template<class LO, class SC>
class BEBilinearFormLaplace1Layer : public BEBilinearForm<LO, SC> {
typedef typename GetType<LO, SC>::SCVT SCVT;
public:
using BEBilinearForm<LO, SC>::assemble;
//! default constructor
BEBilinearFormLaplace1Layer( );
//! copy constructor
BEBilinearFormLaplace1Layer(
const BEBilinearFormLaplace1Layer& orig
);
/*!
* Constructor taking the boundary element space on which the form acts
*/
BEBilinearFormLaplace1Layer(
BESpace<LO, SC>* space,
int* quadratureOrder = nullptr,
quadratureType quadrature = SauterSchwab,
int* quadratureOrderDisjointElems = nullptr,
bool employSymmetricity = false
);
//! destructor
virtual ~BEBilinearFormLaplace1Layer( );
/*!
* The method assembles the Galerkin matrix for the Laplace single layer operator
*/
virtual void assemble(
FullMatrix<LO, SC> & matrix
) const;
/*!
* The method assembles the Galerkin matrix for the Laplace single layer operator
*/
void assembleP0P0(
FullMatrix<LO, SC> & matrix
) const;
void assembleP1P1(
FullMatrix<LO, SC> & matrix
) const;
void assembleP1DisP1Dis(
FullMatrix<LO, SC> & matrix
) const;
protected:
private:
//! default quadrature rule for this class (specific values set in .cpp file)
static int defaultQuadratureOrder[2];
};
}
// include .cpp file to overcome linking problems due to templates
#include "BEBilinearFormLaplace1Layer.cpp"
#endif /* BEBILINEARFORMLAPLACE1LAYER_H */
/*!
* @file BEBilinearFormLaplace2Layer.cpp
* @author Michal Merta
* @date August 8, 2013
*
*/
#ifdef BEBILINEARFORMLAPLACE2LAYER_H
namespace bem4i {
template<class LO, class SC>
BEBilinearFormLaplace2Layer<LO, SC>::BEBilinearFormLaplace2Layer( ) {
}
template<class LO, class SC>
BEBilinearFormLaplace2Layer<LO, SC>::BEBilinearFormLaplace2Layer(
const BEBilinearFormLaplace2Layer& orig
) {
this->space = orig.space;
this->quadratureOrder = orig.quadratureOrder;
this->quadratureOrderDisjointElems = orig.quadratureOrderDisjointElems;
}
template<class LO, class SC>
BEBilinearFormLaplace2Layer<LO, SC>::BEBilinearFormLaplace2Layer(
BESpace<LO, SC>* space,
int* quadratureOrder,
quadratureType quadrature,
int* quadratureOrderDisjointElems
) {
this->quadrature = quadrature;
this->space = space;
if ( quadratureOrder ) {
this->quadratureOrder = quadratureOrder;
} else {
switch ( quadrature ) {
case SauterSchwab:
this->quadratureOrder = defaultQuadraturesSauterSchwab;
break;
case Steinbach:
this->quadratureOrder = defaultQuadraturesSteinbach;
}
}
this->quadratureOrderDisjointElems = quadratureOrderDisjointElems;
}
template<class LO, class SC>
BEBilinearFormLaplace2Layer<LO, SC>::~BEBilinearFormLaplace2Layer( ) {
}
template<class LO, class SC>
void BEBilinearFormLaplace2Layer<LO, SC>::assemble(
FullMatrix<LO, SC>& matrix
) const {
if ( this->space->getAnsatzFunctionType( ) == p1 &&
this->space->getTestFunctionType( ) == p0 ) {
this->assembleP0P1( matrix );
return;
} else if ( this->space->getAnsatzFunctionType( ) == p1 &&
this->space->getTestFunctionType( ) == p1 ) {
this->assembleP1P1( matrix );
return;
} else if ( this->space->getAnsatzFunctionType( ) == p0 &&
this->space->getTestFunctionType( ) == p0 ) {
this->assembleP0P0( matrix );
return;
} else if ( this->space->getTestFunctionType( ) == p1dis &&
this->space->getAnsatzFunctionType( ) == p1dis ) {
this->assembleP1DisP1Dis( matrix );
return;
}
}
template<class LO, class SC>
void BEBilinearFormLaplace2Layer<LO, SC>::assembleP0P1(
FullMatrix<LO, SC> & matrix
) const {
LO nNodes = this->space->getMesh( )->getNNodes( );
LO nElems = this->space->getMesh( )->getNElements( );
matrix.resize( nElems, nNodes, false );
SC * matrixData = matrix.getData( );
#pragma omp parallel for schedule( static, 32 )
for ( LO j = 0; j < nNodes; ++j ) {
for ( LO i = 0; i < nElems; ++i ) {
matrixData[ j * nElems + i ] = 0.0;
}
}
#pragma omp parallel
{
BEIntegratorLaplace<LO, SC> integrator( this->space, this->quadratureOrder,
this->quadrature, this->quadratureOrderDisjointElems );
FullMatrix< LO, SC > elemMatrix( 1, 3 );
LO elem[ 3 ];
#pragma omp for schedule( dynamic, 8 )
for ( LO i = 0; i < nElems; ++i ) {
for ( LO j = 0; j < nElems; ++j ) {
integrator.getElemMatrix2Layer( i, j, elemMatrix );
this->space->getMesh( )->getElement( j, elem );
matrixData[ elem[ 0 ] * nElems + i ] += elemMatrix.get( 0, 0 );
matrixData[ elem[ 1 ] * nElems + i ] += elemMatrix.get( 0, 1 );
matrixData[ elem[ 2 ] * nElems + i ] += elemMatrix.get( 0, 2 );
}
}
} // end omp parallel
}
template<class LO, class SC>
void BEBilinearFormLaplace2Layer<LO, SC>::assembleP1P1(
FullMatrix<LO, SC> & matrix
) const {
LO nNodes = this->space->getMesh( )->getNNodes( );
LO nElems = this->space->getMesh( )->getNElements( );
matrix.resize( nNodes, nNodes, false );
SC * matrixData = matrix.getData( );
#pragma omp parallel for schedule( static, 32 )
for ( LO j = 0; j < nNodes; ++j ) {
for ( LO i = 0; i < nNodes; ++i ) {
matrixData[ j * nNodes + i ] = 0.0;
}
}
#pragma omp parallel
{
BEIntegratorLaplace<LO, SC> integrator( this->space,
this->quadratureOrder, this->quadrature,
this->quadratureOrderDisjointElems );
FullMatrix< LO, SC > elemMatrix( 3, 3 );
LO iElem[ 3 ], oElem[ 3 ];
SC * elemMatrixData = elemMatrix.getData( );
#pragma omp for schedule( dynamic, 32 )
for ( LO i = 0; i < nElems; ++i ) { // i outer row
for ( LO j = 0; j < nElems; ++j ) { // j inner col
integrator.getElemMatrix2Layer( i, j, elemMatrix );
this->space->getMesh( )->getElement( i, oElem );
this->space->getMesh( )->getElement( j, iElem );
for ( int oRot = 0; oRot < 3; ++oRot ) {
for ( int iRot = 0; <