Feel++ 0.91.0
Public Types | Public Member Functions

ResidualEstimator< Dim, Order > Class Template Reference

#include <residualestimator.hpp>

Inheritance diagram for ResidualEstimator< Dim, Order >:
Feel::Simget

List of all members.

Public Types

typedef double value_type
 numerical type is double
typedef Backend< value_typebackend_type
 linear algebra backend factory
typedef boost::shared_ptr
< backend_type
backend_ptrtype
 linear algebra backend factory shared_ptr<> type
typedef
backend_type::sparse_matrix_type 
sparse_matrix_type
 sparse matrix type associated with backend
typedef
backend_type::sparse_matrix_ptrtype 
sparse_matrix_ptrtype
 sparse matrix type associated with backend (shared_ptr<> type)
typedef backend_type::vector_type vector_type
 vector type associated with backend
typedef
backend_type::vector_ptrtype 
vector_ptrtype
 vector type associated with backend (shared_ptr<> type)
typedef Simplex< Dim > convex_type
 geometry entities type composing the mesh, here Simplex in Dimension Dim of Order 1
typedef Mesh< convex_typemesh_type
 mesh type
typedef boost::shared_ptr
< mesh_type
mesh_ptrtype
 mesh shared_ptr<> type
typedef FunctionSpace
< mesh_type, bases< Lagrange
< 0, Scalar, Discontinuous > > > 
p0_space_type
 function space that holds piecewise constant ( $P_0$) functions (e.g. to store material properties or partitioning
typedef boost::shared_ptr
< p0_space_type
p0_space_ptrtype
typedef p0_space_type::element_type p0_element_type
 an element type of the $P_0$ discontinuous function space
typedef bases< Lagrange
< 1, Scalar > > 
p1_basis_type
 the basis type of our approximation space
typedef FunctionSpace
< mesh_type, p1_basis_type
p1_space_type
typedef boost::shared_ptr
< p1_space_type > 
p1_space_ptrtype
typedef p1_space_type::element_type p1_element_type
typedef bases< Lagrange< Order,
Scalar > > 
basis_type
 the basis type of our approximation space
typedef FunctionSpace
< mesh_type, basis_type
space_type
 the approximation function space type
typedef boost::shared_ptr
< space_type
space_ptrtype
 the approximation function space type (shared_ptr<> type)
typedef space_type::element_type element_type
 an element type of the approximation function space
typedef Exporter< mesh_type, 1 > export_type
 the exporter factory type
typedef boost::shared_ptr
< export_type
export_ptrtype
 the exporter factory (shared_ptr<> type)

Public Member Functions

 ResidualEstimator (AboutData const &about)
 ResidualEstimator (po::variables_map const &vm, AboutData const &about)
void run ()
void run (const double *X, unsigned long P, double *Y, unsigned long N)
 BOOST_PARAMETER_CONST_MEMBER_FUNCTION ((mesh_ptrtype), adapt, tag,(required(h,*))(optional(maxit,*(boost::is_integral< mpl::_ >), 10)(hmin,*(boost::is_arithmetic< mpl::_ >), 1e-2)(hmax,*(boost::is_arithmetic< mpl::_ >), 2)(model,*,"")(statistics,*(boost::is_integral< mpl::_ >), 0)(update,*(boost::is_integral< mpl::_ >), MESH_CHECK|MESH_UPDATE_FACES|MESH_UPDATE_EDGES|MESH_RENUMBER)(collapseOnBoundary,*(boost::is_integral< mpl::_ >), true)(collapseOnBoundaryTolerance,*(boost::is_arithmetic< mpl::_ >), 1e-6)))

Detailed Description

template<int Dim, int Order = 1>
class ResidualEstimator< Dim, Order >

Laplacian Solver using continuous approximation spaces solve $ -\Delta u = f$ on $\Omega$ and $u= g$ on $\Gamma$

Template Parameters:
Dimthe geometric dimension of the problem (e.g. Dim=1, 2 or 3)
Orderthe approximation order

Constructor & Destructor Documentation

template<int Dim, int Order = 1>
ResidualEstimator< Dim, Order >::ResidualEstimator ( AboutData const &  about) [inline]

Constructor


Member Function Documentation

template<int Dim, int Order>
void ResidualEstimator< Dim, Order >::run ( ) [virtual]

simply execute the simget

Implements Feel::Simget.

template<int Dim, int Order>
void ResidualEstimator< Dim, Order >::run ( const double *  X,
unsigned long  P,
double *  Y,
unsigned long  N 
) [virtual]

models the input/output relation $ Y=F(X) $

The function space and some associated elements(functions) are then defined

 */
    P0h = p0_space_type::New( mesh );
    P1h = p1_space_type::New( mesh );
    space_ptrtype Xh = space_type::New( mesh );
    element_type u( Xh, "u" );
    element_type v( Xh, "v" );

define $g$ the expression of the exact solution and $f$ the expression of the right hand side such that $g$ is the exact solution

 */
    //# marker1 #
    value_type pi = M_PI;
    auto fn1 = (1-Px()*Px())*(1-Py()*Py())*(1-Pz()*Pz())*exp(beta*Px());
    auto fn2 = sin(alpha*pi*Px())*cos(alpha*pi*Py())*cos(alpha*pi*Pz())*exp(beta*Px());
    auto g = chi(fn==1)*fn1 + chi(fn==2)*fn2;
    auto grad_g =
        trans(chi(fn==1)*((-2*Px()*(1-Py()*Py())*(1-Pz()*Pz())*exp(beta*Px())+beta*fn1)*unitX()+
                          -2*Py()*(1-Px()*Px())*(1-Pz()*Pz())*exp(beta*Px())*unitY()+
                          -2*Pz()*(1-Px()*Px())*(1-Py()*Py())*exp(beta*Px())*unitZ())+
              chi(fn==2)*(+(alpha*pi*cos(alpha*pi*Px())*cos(alpha*pi*Py())*cos(alpha*pi*Pz())*exp(beta*Px())+beta*fn2)*unitX()+
                          -alpha*pi*sin(alpha*pi*Px())*sin(alpha*pi*Py())*cos(alpha*pi*Pz())*exp(beta*Px())*unitY()+
                          -alpha*pi*sin(alpha*pi*Px())*cos(alpha*pi*Py())*sin(alpha*pi*Pz())*exp(beta*Px())*unitZ()));
    auto minus_laplacian_g =
        (chi( fn == 1 )*( 2*(1-Py()*Py())*(1-Pz()*Pz())*exp(beta*Px()) + 4*beta*Px()*(1-Py()*Py())*(1-Pz()*Pz())*exp(beta*Px()) - beta*beta *fn1 +
                          2*(1-Px()*Px())*(1-Pz()*Pz())*exp(beta*Px())*chi(Dim >= 2) +
                          2*(1-Px()*Px())*(1-Py()*Py())*exp(beta*Px())*chi(Dim == 3) )  +
         chi( fn == 2 )*  (
             exp(beta*Px())*(Dim*alpha*alpha*pi*pi*sin(alpha*pi*Px())-beta*beta*sin(alpha*pi*Px())-2*beta*alpha*pi*cos(alpha*pi*Px()))*
             ( cos(alpha*pi*Py())*chi(Dim>=2) + chi(Dim==1)) * ( cos(alpha*pi*Pz())*chi(Dim==3) + chi(Dim<=2) )
             )

            );

    //# endmarker1 #

Construction of the right hand side. F is the vector that holds the algebraic representation of the right habd side of the problem

 */
    //# marker2 #
    vector_ptrtype F( M_backend->newVector( Xh ) );
    form1( _test=Xh, _vector=F, _init=true ) =
        integrate( elements(mesh), minus_laplacian_g*id(v) )+
        integrate( markedfaces( mesh, tag_Neumann ),
                   grad_g*vf::N()*id(v) );
    //# endmarker2 #
    if ( this->comm().size() != 1 || weakdir )
    {
        //# marker41 #
        form1( _test=Xh, _vector=F ) +=
            integrate( markedfaces(mesh,tag_Dirichlet), g*(-grad(v)*vf::N()+penaldir*id(v)/hFace()) );
        //# endmarker41 #
    }
    F->close();

create the matrix that will hold the algebraic representation of the left hand side

 */
    sparse_matrix_ptrtype D( M_backend->newMatrix( Xh, Xh ) );

assemble $ u v$

 */
    form2( Xh, Xh, D, _init=true ) =
        integrate( elements(mesh), gradt(u)*trans(grad(v)) );

weak dirichlet conditions treatment for the boundaries marked 1 and 3

  1. assemble $\int_{\partial \Omega} -\nabla u \cdot \mathbf{n} v$
  2. assemble $\int_{\partial \Omega} -\nabla v \cdot \mathbf{n} u$
  3. assemble $\int_{\partial \Omega} \frac{\gamma}{h} u v$
 */
        //# marker10 #
        form2( Xh, Xh, D ) +=
            integrate( markedfaces(mesh,tag_Dirichlet),
                       -(gradt(u)*vf::N())*id(v)
                       -(grad(v)*vf::N())*idt(u)
                       +penaldir*id(v)*idt(u)/hFace());
        D->close();
        //# endmarker10 #

strong(algebraic) dirichlet conditions treatment for the boundaries marked 1 and 3

  1. first close the matrix (the matrix must be closed first before any manipulation )
  2. modify the matrix by cancelling out the rows and columns of D that are associated with the Dirichlet dof
 */
        //# marker5 #
        D->close();
        form2( Xh, Xh, D ) +=
            on( markedfaces(mesh, tag_Dirichlet), u, F, g );
        //# endmarker5 #

solve the system

 */
    backend_type::build()->solve( _matrix=D, _solution=u, _rhs=F );

*/

//! compute the

Implements Feel::Simget.

References Feel::Environment::changeRepository(), Feel::element_product(), Feel::elements(), Feel::internalfaces(), and Feel::markedfaces().