Posted by
Nicolas Di Césaré on
Oct 05, 2001; 7:28am
URL: http://quantlib.414.s1.nabble.com/Docs-requesting-peer-review-tp1768p1771.html
En réponse à Luigi Ballabio <
[hidden email]>:
>
> The work is done jointly by three methods, namely,
> TridiagonalOperatorCommon::setLowerBC(bc)/setUpperBC(bc) and
> TridiagonalOperatorCommon::solveFor(rhs).
>
> The former store the value of the boundary condition and, in the
> Dirichlet
> case, set the upper row of the operator to [1,0,0,....] or the lower row
> to
> [0,0,....,0,1], respectively.
>
> The latter sets the first/last element of the rhs to the stored boundary
>
> condition value(s), then solves the tridiagonal system. This, together
> with
> the operator modification made previously, ensures that the result
> satisfies the boundary condition.
>
> Similar tricks are done in the Neumann case.
>
> The drawback I see for the current implementation is that if one creates
> an
> operator with a built-in boundary condition and then sets a Dirichlet or
>
> Neumann one, then it is not possible to revert to the original b.c.
> unless
> one manually sets the first/last row of the operator. With a bit of
> housekeeping, this could be done in a future implementation. However,
> I'm
> not at all certain whether the reverting operation should be triggered
> by
> setting a b.c. with type BoundaryCondition::None (meaning "do nothing to
> the operator I passed"), or rather some other type we should add to the
> enumeration (meaning "the original b.c. the operator had" as opposed to
> None meaning "no b.c.").
> Feel free to make any suggestions---or even to implement it :)
I have worked on a similar design for finite element
method applied to stationnary equations (like Laplace or Poisson equations)
that can be view as backward Euler equation. My conclusion was that boundary
conditions (BC) can modify the equation operator and the right hand side (RHS)
of the equation. Thus, the BC base class interface is :
template <class Vec, class Mat>
class BaseBdyCondition {
public:
virtual
void applyTo(const Vec& bdyval, // BC values
Mat & sparsemat, // Matrix representation of the operator
Vec & systemrhs) // RHS
= 0;
};
With that interface, I have implemented Neumann, Dirichlet(penalization method)
and Fourrier BC.
The resolution of the problem is made as follow :
// Define problem type
Elliptic_PDE Problem(Mesh);
// Build the equation operator ie set the matrix
Problem.Equation(laplacian-id);
// Set the right hand side
Problem.Set_Rhs(Source);
// Set Bdy Cond type
Problem.Bdy_Conditions_Types(Flow_Bdy_Cond);
// Set Bdy Cond value
Problem.Bdy_Conditions(Dirichlet_Bdy_Values);
// Set the resolution method
Problem.Resolution_Method(gauss);
// solve the linear system
Problem.Solve();
This method is wrong for explicit problem where the BC is treated after the
time step u^(k-1) = A.u^k
In conclusion, if we want to allow generic BC to be implemented we have to
1 - Separate resolution of explicit and implicit schemes.
2 - use a more generic BC design (like mine ;-) )
--
Nicolas Di Césaré
http://acm.emath.fr/~dicesare