packages/Alea/mc_solvers/LinearSystemFactory.cc
//---------------------------------------------------------------------------//
// Build cell-based 1D finite difference Laplacian matrix on [0,1]
//---------------------------------------------------------------------------//
Build cell-based 1D finite difference Laplacian matrix on [0,1]
Build cell - based 1D finite difference Laplacian matrix on [ 0 1 ]
Teuchos::RCP<CRS_MATRIX>
LinearSystemFactory::buildLaplacianMatrix(
    int N, Teuchos::RCP<Teuchos::ParameterList> pl )
{
    SCALAR h = 1.0 / static_cast<SCALAR>(N);
    SCALAR laplacian_shift = pl->get("laplacian_shift",0.0);
    VALIDATE( laplacian_shift >= 0.0, "laplacian shift must be non-negative");
    Teuchos::RCP<const Teuchos::Comm<int> > comm =
        Teuchos::DefaultComm<int>::getComm();
    LO base=0;
    Teuchos::RCP<MAP> map( new MAP(N,base,comm) );
    LO numPerRow = 3;
    Teuchos::RCP<CRS_MATRIX> A =
        Teuchos::rcp( new CRS_MATRIX( map, numPerRow ) );
    LO numMyRows = map->getNodeNumElements();
    Teuchos::ArrayRCP<SCALAR> val(3);
    Teuchos::ArrayRCP<GO>     ind(3);
    for( LO i=0; i<numMyRows; ++i )
    {
        LO numThisRow;
        GO gid = map->getGlobalElement( i );
        if( gid == 0 )
        {
            ind[0] = 0;
            val[0] = (3.0+laplacian_shift)/(h*h);
            ind[1] = 1;
            val[1] = -1.0/(h*h);
            numThisRow = 2;
        }
        else if( gid == N-1 )
        {
            ind[0] = N-2;
            val[0] = -1.0/(h*h);
            ind[1] = N-1;
            val[1] = (3.0+laplacian_shift)/(h*h);
            numThisRow = 2;
        }
        else
        {
            ind[0] = gid-1;
            val[0] = -1.0/(h*h);
            ind[1] = gid;
            val[1] = (2.0+laplacian_shift)/(h*h);
            ind[2] = gid+1;
            val[2] = -1.0/(h*h);
            numThisRow = 3;
        }
        Teuchos::ArrayView<SCALAR> val_view = val(0,numThisRow);
        Teuchos::ArrayView<GO>     ind_view = ind(0,numThisRow);
        A->insertGlobalValues( gid, ind_view, val_view );
    }
    A->fillComplete();
    return A;
}
Teuchos :: RCP < CRS_MATRIX > LinearSystemFactory :: buildLaplacianMatrix ( int N , Teuchos :: RCP < Teuchos :: ParameterList > pl ) { SCALAR h = 1.0 / static_cast < SCALAR > ( N ) ; SCALAR laplacian_shift = pl -> get ( " " , 0.0 ) ; VALIDATE ( laplacian_shift >= 0.0 , " " ) ; Teuchos :: RCP < const Teuchos :: Comm < int > > comm = Teuchos :: DefaultComm < int > :: getComm ( ) ; LO base = 0 ; Teuchos :: RCP < MAP > map ( new MAP ( N , base , comm ) ) ; LO numPerRow = 3 ; Teuchos :: RCP < CRS_MATRIX > A = Teuchos :: rcp ( new CRS_MATRIX ( map , numPerRow ) ) ; LO numMyRows = map -> getNodeNumElements ( ) ; Teuchos :: ArrayRCP < SCALAR > val ( 3 ) ; Teuchos :: ArrayRCP < GO > ind ( 3 ) ; for ( LO i = 0 ; i < numMyRows ; ++ i ) { LO numThisRow ; GO gid = map -> getGlobalElement ( i ) ; if ( gid == 0 ) { ind [ 0 ] = 0 ; val [ 0 ] = ( 3.0 + laplacian_shift ) / ( h * h ) ; ind [ 1 ] = 1 ; val [ 1 ] = -1.0 / ( h * h ) ; numThisRow = 2 ; } else if ( gid == N - 1 ) { ind [ 0 ] = N - 2 ; val [ 0 ] = -1.0 / ( h * h ) ; ind [ 1 ] = N - 1 ; val [ 1 ] = ( 3.0 + laplacian_shift ) / ( h * h ) ; numThisRow = 2 ; } else { ind [ 0 ] = gid - 1 ; val [ 0 ] = -1.0 / ( h * h ) ; ind [ 1 ] = gid ; val [ 1 ] = ( 2.0 + laplacian_shift ) / ( h * h ) ; ind [ 2 ] = gid + 1 ; val [ 2 ] = -1.0 / ( h * h ) ; numThisRow = 3 ; } Teuchos :: ArrayView < SCALAR > val_view = val ( 0 , numThisRow ) ; Teuchos :: ArrayView < GO > ind_view = ind ( 0 , numThisRow ) ; A -> insertGlobalValues ( gid , ind_view , val_view ) ; } A -> fillComplete ( ) ; return A ; }
Build cell-based 1D finite difference Laplacian matrix on [0,1]
Build cell - based 1D finite difference Laplacian matrix on [ 0 1 ]
// Build comm and map // Dirichlet on left // Dirichlet on right