Commit 08e7ba7c authored by Jan Zapletal's avatar Jan Zapletal

added regularization vector

parent 7985d943
......@@ -47,7 +47,9 @@ int main(
void testHeatDtN(
) {
double * Sarray = new double [ nNodes * nNodes ];
double * S = new double [ nNodes * nNodes ];
double * a = new double [ nNodes ];
double beta;
double alpha = 1.0;
int qType = 1;
int orderNear = 4;
......@@ -58,7 +60,9 @@ void testHeatDtN(
bool nonsymmetric = false;
bool verbose = true;
bem4i::getLaplaceSteklovPoincare( Sarray,
bem4i::getLaplaceSteklovPoincare( S,
a,
&beta,
nNodes,
nodes,
nElems,
......@@ -73,15 +77,28 @@ void testHeatDtN(
bem4i::deleteBem4iData<int, double>( data );
if( verbose )
if( verbose ){
std::cout << "S: " << std::endl;
for( size_t i = 0; i < nNodes; ++i ){
for( size_t j = 0; j < nNodes; ++j ){
std::cout << Sarray[ j * nNodes + i ] << " ";
std::cout << S[ j * nNodes + i ] << " ";
}
std::cout << std::endl;
}
std::cout << "a: " << std::endl;
for( size_t i = 0; i < nNodes; ++i ){
std::cout << a[ i ] << " ";
}
std::cout << std::endl;
std::cout << "beta: " << std::endl;
std::cout << beta << std::endl;
}
delete [] Sarray;
delete [] S;
delete [] a;
}
......@@ -52,7 +52,9 @@ void deleteBem4iData(
template<class LO, class SC>
void getLaplaceSteklovPoincare(
SC * Sarray,
SC * S,
SC * a,
SC * beta,
LO nNodes,
const SC * nodes,
LO nElems,
......@@ -132,6 +134,21 @@ void getLaplaceSteklovPoincare(
SparseMatrix< LO, SC > M;
id.assemble( M );
if( a != nullptr ){
Vector< LO, SC > avect( nNodes, a, false );
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;
}
}
Eigen::SparseMatrix<SC, Eigen::ColMajor, LO> * Me = M.getEigenSparseMatrix( );
// K <- 1/2M + K
......@@ -143,18 +160,18 @@ void getLaplaceSteklovPoincare(
}
}
FullMatrix< LO, SC > S( nNodes, nNodes, Sarray );
FullMatrix< LO, SC > Smat( nNodes, nNodes, S );
BEBilinearFormLaplaceHypersingular< LO, SC > formD( &bespace11, quadNear,
quadType, quadFar );
if ( !nonsymmetric ) {
// S <- D
if ( verbose ) ProgressMonitor::init( "D" );
formD.assemble( S, *V );
//formD.assemble( S );
formD.assemble( Smat, *V );
//formD.assemble( Smat );
if ( verbose ) ProgressMonitor::step( );
} else {
// S <- 0
S.setAll( 0.0 );
Smat.setAll( 0.0 );
}
FullMatrix< LO, SC > * VinvK = new FullMatrix< LO, SC >( *K );
......@@ -165,13 +182,13 @@ void getLaplaceSteklovPoincare(
if ( !nonsymmetric ) {
// S <- D + (1/2M + K)' * V^-1 * (1/2M + K)
S.multiply( *K, *VinvK, true, false, 1.0, 1.0 );
Smat.multiply( *K, *VinvK, true, false, 1.0, 1.0 );
} else {
// S <- M' * V^-1 * (1/2M + K)
S.multiply( M, *VinvK, true, false, 1.0, 0.0 );
Smat.multiply( M, *VinvK, true, false, 1.0, 0.0 );
}
S.scale(alpha);
Smat.scale(alpha);
delete VinvK;
data->V = V;
......@@ -240,7 +257,9 @@ template void bem4i::deleteBem4iData< int, double >(
);
template void bem4i::getLaplaceSteklovPoincare< int, double >(
double * Sarray,
double * S,
double * a,
double * beta,
int nNodes,
const double * nodes,
int nElems,
......@@ -277,7 +296,9 @@ template void bem4i::deleteBem4iData< long, double >(
);
template void bem4i::getLaplaceSteklovPoincare< long, double >(
double * Sarray,
double * S,
double * a,
double * beta,
long nNodes,
const double * nodes,
long nElems,
......
......@@ -12,7 +12,9 @@ void deleteBem4iData(
template<class LO, class SC>
void getLaplaceSteklovPoincare(
SC * Sarray,
SC * S,
SC * a,
SC * beta,
LO nNodes,
const SC * nodes,
LO nElems,
......
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