run the convergence test
The function space and some associated elements(functions) are then defined
*/
space_ptrtype Xh = space_type::New( mesh );
Xh->printInfo();
element_type u( Xh, "u" );
element_type v( Xh, "v" );
element_type gproj( 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
*/
bool weak_dirichlet = this->vm()["weakdir"].template as<int>();
value_type penaldir = this->vm()["penaldir"].template as<double>();
value_type nu = this->vm()["nu"].template as<double>();
std::string exact =
this->vm()[(boost::format("exact%1%D")%Dim).str()].template
as<std::string>();
std::string rhs =
this->vm()[(boost::format("rhs%1%D")%Dim).str()].template
as<std::string>();
LOG(INFO) << "exact = " << exact << "\n";
auto vars=symbols<Dim>();
ex ff;
if ( rhs.empty() )
ff=-nu*laplacian(exact,vars);
else
ff = parse(rhs,vars);
LOG(INFO) << "rhs="<< ff << "\n";
auto g = expr(exact,vars);
auto f = expr(ff,vars);
auto gradg = expr<1,Dim,2>(grad(exact,vars), vars );
LOG(INFO) << "grad(g)="<< grad(exact,vars) << "\n";
gproj = project( _space=Xh, _range=elements( mesh ), _expr=g );
Construction of the right hand side. F is the vector that holds the algebraic representation of the right habd side of the problem
*/
auto F = backend()->newVector( Xh );
auto l = form1( _test=Xh, _vector=F );
l = integrate( _range=elements( mesh ), _expr=f*id( v ) )+
integrate( _range=markedfaces( mesh, "Neumann" ),
_expr=nu*gradg*vf::N()*id( v )
);
if ( weak_dirichlet )
{
l += integrate( _range=markedfaces( mesh,"Dirichlet" ),
_expr=nu*g*( -grad( v )*vf::N()
+ penaldir*id( v )/hFace() ) );
}
create the matrix that will hold the algebraic representation of the left hand side
*/
auto D = backend()->newMatrix( _test=Xh, _trial=Xh );
assemble $ u v$
*/
auto a = form2( _test=Xh, _trial=Xh, _matrix=D );
a = integrate( _range=elements( mesh ), _expr=nu*gradt( u )*trans( grad( v ) ) );
weak dirichlet conditions treatment for the boundaries marked 1 and 3
- assemble \(\int_{\partial \Omega} -\nabla u \cdot \mathbf{n} v\)
- assemble \(\int_{\partial \Omega} -\nabla v \cdot \mathbf{n} u\)
- assemble \(\int_{\partial \Omega} \frac{\gamma}{h} u v\)
*/
a += integrate( _range=markedfaces( mesh,"Dirichlet" ),
_expr= nu * ( -( gradt( u )*vf::N() )*id( v )
-( grad( v )*vf::N() )*idt( u )
+penaldir*id( v )*idt( u )/hFace() ) );
strong(algebraic) dirichlet conditions treatment for the boundaries marked 1 and 3
- first close the matrix (the matrix must be closed first before any manipulation )
- modify the matrix by cancelling out the rows and columns of D that are associated with the Dirichlet dof
*/
a += on( _range=markedfaces( mesh, "Dirichlet" ),
_element=u, _rhs=F, _expr=g );
solve the system
*/
backend( _rebuild=true )->solve( _matrix=D, _solution=u, _rhs=F );
compute the
*/
double L2error =normL2( _range=elements( mesh ),_expr=( idv( u )-g ) );
LOG(INFO) << "||error||_L2=" << L2error << "\n";
save the results
*/
element_type e( Xh, "e" );
e = project( _space=Xh, _range=elements( mesh ), _expr=g );
export_ptrtype exporter( export_type::New() );
if ( exporter->doExport() )
{
LOG(INFO) << "exportResults starts\n";
exporter->step( 0 )->setMesh( mesh );
exporter->step( 0 )->add( "solution", u );
exporter->step( 0 )->add( "exact", e );
exporter->save();
LOG(INFO) << "exportResults done\n";
}